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EXECUTIVE  SUMMARY 

This  document  presents  the  results  of  a  24  month  contract  for  data- 
intensive  analysis  of  the  amplitudes,  travel  times,  and  coda  characteristics 
of  teleseismic  short  period  P  waves  from  underground  explosions.  The  basic 
objective  of  this  line  of  research  is  to  attain  an  understanding  of  the 
near-source  contributions  to  the  variability  of  the  seismic  observations. 

This  will  improve  our  fundamental  understanding  of  the  seismic  radiation 
and  provide  guidance  for  improved  yield  estimation  methodologies.  In 
addition,  analysis  of  the  scattered  wavefield  in  the  P  coda  will  lead  to 
a  better  knowledge  of  scattering  processes  which  dominate  regional  seismograms. 
The  most  basic  contribution  of  this  work  has  been  the  demonstration  that 
even  using  teleseismic  signals  alone,  it  is  possible  to  characterize  near 
source  effects  on  amplitudes  and  travel  times  for  events  in  a  given  test 
site.  The  importance  of  upper  mantle  velocity  heterogeneity  is  strongly 
demonstrated,  as  is  the  relative  stability  of  P  wave  coda.  We  have  made 
the  first  attempt  to  deterministically  analyze  near-source  scattering  in 
the  P  wave  coda,  as  well  as  statistically  documenting  magnitude  dependent 
relative  anomalies  between  direct  arrivals  and  P  coda. 
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Modeling  Teleseismic  SV  Waves  From  Underground  Explosions 
With  Tectonic  Release:  Results  for  Southern  Novaya  Zemlya 
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MODELING  TELESEISMIC  SV  WAVES  FROM  UNDERGROUND 
EXPLOSIONS  WITH  TECTONIC  RELEASE:  RESULTS  FOR  SOUTHERN 

NOVAYA  ZEMLYA 

By  Brian  P.  Cohee  and  Thorne  Lay 
Abstract 

Detailed  forward  modeling  of  long-period  shear  waves  for  two  large  un¬ 
derground  explosions  at  the  Southern  Novaya  Zemlya  test  site  indicates  that 
the  appropriate  equivalent  doubie-couple  orientation  for  the  tectonic  release 
radiation  is  vertical  strike-slip.  Previous  studies  of  observed  teleseismic  SH 
waveforms  and  SV  amplitudes  for  the  27  October  1973  and  2  November  1974 
events  using  geometric  ray  theory  could  not  distinguish  between  vertical 
strike-slip  and  45°-dipping  thrust  geometries.  Either  mechanism  can  match 
the  observed  four-lobed  SH  radiation  pattern,  and  the  two-lobed  SV  ampli¬ 
tude  pattern  can  be  produced  by  interference  with  an  appropriate  size  explo¬ 
sion  pS  signal.  However,  the  complexity  of  the  observed  SV  waveforms  arising 
from  Sp  conversions  near  the  receiver,  diffracted  Sp,  and  shear-coupled  PL 
phases  is  no!  accounted  for  in  the  ray  theory  synthetics.  Incorporating  more 
realistic  Green's  functions  using  Baag  and  Langston's  (1985b)  WKBJ  spectral 
method  allows  more  complete  modeling  of  the  SV  signals.  Due  to  differences 
in  frequency  content  between  the  explosion  and  double-couple  SV  waveforms, 
constructive  interference  occurs  more  efficiently  than  destructive  interfer¬ 
ence  when  the  two  signals  are  linearly  superimposed.  As  a  result,  using  tec¬ 
tonic  release  moments  determined  from  the  SH  waves  and  the  optimum  F 
factors  required  to  match  the  SV  amplitude  patterns,  the  waveforms  produced 
by  the  strike-slip  and  thrust  orientations  differ  substantially  at  some  azimuths. 

The  strike-slip  solution  yields  a  consistently  superior  match  to  the  data.  Using 
the  EU2  model  of  Lerner-Lam  and  Jordan  (1987)  for  the  source  region  and 
either  EU2  or  TNA  (Grand  and  Helmberger,  1984)  for  the  receiver  structure, 
together  with  an  attenuation  model  similar  to  SL8,  we  obtain  a  doubie-couple 
moment,  M0  =  3.2  x  1024  dyne-cm  and  explosion  source  strength,  Jj*  =  3.8 
±  0.5  x  10"  cm3  for  the  27  October  1973  event,  and  M0  =  1.7  x  1024  dyne- 
cm  and  =  2.0  ±  0.3  x  10"  cm3  for  the  2  November  1974  event.  Complete 
waveform  modeling  of  SV  signals  can  thus  provide  improved  constraints  on 
tectonic  release  radiation  and  explosion  source  strength. 

Introduction 

Numerous  studies  have  addressed  the  existence  of  nonisotropic  seismic  radia¬ 
tion  accompanying  large  underground  nuclear  explosions  (e.g.  Press  and  Ar- 
chambeau,  1962;  Toksoz  et  al.,  1965;  Hirasawa,  1971;  Wallace  et  al.,  1983,  1985). 
Long-period  horizontal  shear  wave  motions,  including  Love  waves  and  tele¬ 
seismic  SH  body  waves,  result  directly  from  this  nonisotropic  component  and  can 
be  used  to  constrain  the  transverse  source  contribution.  Although  many  uncer¬ 
tainties  remain  concerning  the  precise  physical  source  of  this  transverse  energy, 
it  is  generally  attributed  to  the  release  of  preexisting  tectonic  strain  accumulation 
in  the  vicinity  of  the  explosion  (tectonic  release).  Two  general  theories  for  tectonic 
release  have  been  suggested:  the  triggering  of  displacement  on  a  nearby  fault 
IBrune  and  Pomeroy,  1963;  Aki  et  al.,  1969;  Aki  and  Tsai,  1972),  or  the  release 
of  in  situ  stress  from  vaporized  and  fractured  rock  during  formation  of  the  ex- 
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plosion  cavity  (Press  and  Archambeau,  1962;  Archambeau  and  Sammis,  1970; 
Archambeau,  1972;  Lambert  et  al.,  1972).  For  either  tectonic  release  model,  the 
long-period  far-field  radiation  may  be  represented  by  an  equivalent  double-couple 
dislocation;  however,  it  is  challenging  to  find  a  unique  focal  mechanism  solution. 

The  ambiguity  intrinsic  to  determination  of  the  tectonic  release  focal  mecha¬ 
nism  using  long-period  S  waves  was  explored  by  Wallace  et  al.  (1985).  Both  ver¬ 
tical  strike-slip  and  45°-dipping  (thrust  or  normal)  faults  can  produce  identical 
four-lobed  SH  (and  Love  wave)  radiation  patterns,  with  the  corresponding  mo¬ 
ments  for  the  dip-slip  mechanisms  being  doubled  and  strike  azimuths  being  ro¬ 
tated  ±45°  relative  to  the  strike-slip  mechanism  (see  Aki  and  Tsai,  1972;  Given 
and  Mellman,  1385;  Burger  et  al.,  1986b).  As  a  result,  uncertainty  in  the  true 
orientation  has  direct  consequences  for  the  appraisal  of  long-period  source  am¬ 
plitude  and  corresponding  explosion  strength.  Wallace  et  al.  (1985)  demonstrate 
that  this  uncertainty  cannot  be  resolved  using  SH  data  alone,  but  the  inclusion 
of  SV  data  can  yield  a  unique  model.  Tectonic  release  studies  incorporating  both 
SH  and  SV  waves  have  still  been  unable  to  uniquely  resolve  the  double-couple 
orientation  in  cases  (which  are  frequently  encountered)  where  the  focal  mech¬ 
anism  is  either  vertical  strike-slip  or  45°  dip-slip,  because  it  has  been  difficult  to 
quantitatively  model  the  complex  SV  arrivals  (Burger  et  al.,  1986b). 

The  complexity  inherent  in  modeling  SV  wave  propagation  results  from  the 
contributions  of  additional  phases  within  a  30-sec  window  containing  the  direct 
S  arrival.  Sp  conversions  occurring  at  the  Moho  and  other  discontinuities  beneath 
the  receiver  immediately  precede  and  obscure  direct  S.  The  long-period  shear- 
coupled  PL  ( SPL )  phase,  initially  recognized  and  investigated  by  Oliver  (1961), 
follows  and  often  overshadows,  the  initial  pulse  (Baag  and  Langston,  1985a).  The 
SPL  phase  is  set  up  as  a  normally  dispersed  leaking  mode  in  the  crustal  wave 
guide  with  trapped  P  waves  (leaking  downgoing  S  energy)  and  SV  reverberations. 
Diffracted  Sp  (Frazer,  1977)  is  a  diffracted  head  wave  set  up  at  the  base  of  the 
crust  by  a  critically  incident  S  wave,  which  arrives  after  direct  S  at  distances 
greater  than  50°  (Baag  and  Langston,  1986).  The  propagation  paths  and  char¬ 
acteristics  of  these  auxiliary  SV  energy  contributions  are  examined  in  more  detail 
by  Burdick  and  Langston  (1977;  and  Baag  and  Langston  (1985a.  1986'.  These 
additional  phases  are  very  sensitive  to  source  and  receiver  structure,  particularly 
for  shallow  sources,  further  complicating  the  analysis  of  explosion  signals.  The 
difficulty  of  interpreting  the  SV  wave  train  has  limited  its  use  for  source  inves¬ 
tigations.  However,  this  very  complexity  contains  potential  for  increased  reso¬ 
lution  of  the  tectonic  release  orientation  if  the  wave  field  can  be  accurately- 
computed . 

In  forward  modeling  of  teleseismic  signals  from  explosions  with  tectonic  re¬ 
lease,  it  is  commonly  assumed  that  both  sources  are  coincident  in  space  and  time 
and  can  thus  be  represented  by  the  linear  superposition  of  the  explosion  and 
double-couple  waveforms.  This  superposition  necessitates  dependable  Green's 
functions  to  realistically  account  for  interference  between  the  two  source  signals. 

Ray  theory  can  be  used  to  construct  teleseismic  SH  and  SV  synthetics  (Langs¬ 
ton  and  Helmberger,  1975),  which  is  usually  sufficient  to  match  the  simple  SH 
waveforms.  However,  the  SV  synthetics  include  only  the  direct  geometric  phases 
of  the  explosion  ( pS )  and  tectonic  release  (S,  sS,  pSl  signals.  At  distances  less 
than  60°,  these  are  a  poor  approximation  of  the  total  SV  arrival  (Langston  and 
Baag,  1985).  A  more  complete  calculation  that  includes  SPL,  Sp  conversions,  and 
diffracted  Sp,  and  accounts  for  their  excitation  dependence  on  source  type,  is 
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required  to  model  the  complex  interference  between  the  two  sources.  A  significant 
contribution  to  the  understanding  of  the  S  V  system  at  distances  between  30°  and 
90°  was  recently  made  by  Baag  and  Langston  (1985a,  b).  A  WKBJ  spectral 
method  for  creating  teleseismic  synthetic  SV  waves  was  presented,  which  re¬ 
quires  substantially  less  computation  time  than  previously  necessary  to  compute 
adequate  Green’s  functions.  The  WKBJ  spectral  method  couples  WKBJ  theory 
for  a  generalized  ray  in  a  vertically  inhomogeneous  half-space  with  the  Watson 
( 1970)  compound-matrix  technique  in  near-surface  homogeneous  layers  at  both 
the  source  and  receiver.  The  validity  of  the  method  has  been  demonstrated 
through  comparison  with  Kind’s  (1978)  reflectivity  treatment  of  the  P-SV  wave 
propagation  equations  (Baag  and  Langston,  1985b). 

In  this  study,  we  construct  synthetic  SV  waveforms  to  resolve  the  ambiguity 
of  the  double-couple  orientation  of  tectonic  release  for  two  Soviet  underground 
nuclear  tests,  27  October  1973  and  2  November  1974,  and  to  determine  appro¬ 
priate  explosion  and  tectonic  release  source  strengths.  These  two  explosions  are 
among  a  series  of  large  underground  nuclear  tests  conducted  by  the  Soviet  Union 
on  the  island  ofNovaya  Zemlya  in  the  Barents  Sea  (Figure  1).  The  tectonic  release 
focal  mechanisms  for  nine  events  at  a  northern  site  were  found  to  be  oblique 


Fic.  1.  Base  map  showing  the  two  Novava  Zemlva  test  sites.  The  northern  and  southern  test  sites 
are  separated  by  about  300  1cm  The  two  nuclear  explosions  analyzed  in  this  study  occurred  at  the 
southern  site  The  lower  hemisphere  projection  shows  P- wave  radiation  nodes  for  the  two  possible 
tectonic  release  orientations  that  satisfy  the  SH  observations 


6 


MODELING  .SV  WAVES  FROM  EXPLOSIONS  1 161 

normal  by  Burger  et  al.  (1986b)  based  on  the  SH  radiation  pattern.  SH  obser¬ 
vations  from  the  Southern  Novaya  Zemlya  site  however,  could  not  distinguish 
between  vertical  strike-slip  and  45°  dip-slip  double-couple  orientations.  Long- 
period.  teleseismic  SV  waves  for  the  27  October  1973  and  2  November  1974 
Southern  Novaya  Zemlya  events  are  modeled  to  discern  characteristic  traits  of 
each  tectonic  release  double-couple  and  its  interference  with  the  explosion  signal, 
and  to  resolve  the  true  orientation  by  comparison  with  observed  waveforms. 

Shear  Wave  Observations  and  Ray  Theory  Modeling 

The  27  October  1973  and  2  November  1974  Southern  Novaya  Zemlya  explosions 
are  among  the  largest  underground  nuclear  tests  ever  conducted  (see  Table  1). 
As  a  result,  numerous  long-period  WWSSN  and  CSN  records  are  available  from 
teleseismic  distances.  Selected  SH  data  for  the  27  October  1973  event  are  shown 
in  Figure  2.  Systematic  polarity  reversals  near  60'.  150°,  240°,  and  330°  define 
a  four-lobed  amplitude  pattern  with  vertical  nodal  planes.  Both  vertical  strike- 
slip  and  45°  dip-slip  mechanisms  are  characterized  by  such  an  SH  radiation  pat¬ 
tern.  The  2  November  1974  SH  data  are  presented  in  a  similar  format  by  Burger 
et  al.  (1986b,  Figure  3  c  Comparison  of  the  observations  for  both  explosions  reveal 
a  similar  SH  radiation,  indicative  of  a  common  tectonic  release  orientation. 

The  simple  SH  pulses  are  well-modeled  using  ray  theory  for  a  half-space 
(Langston  and  Helmberger,  1975 1  and  a  point-source  dislocation  with  a  3.0-sec 
duration  trapezoid  source-time  function  1 1  sec  rise.  1  sec  top,  1  sec  fall)  convolved 
with  an  appropriate  15-100  WWSSN  instrument  response.  Source  depth  was 
placed  at  1  km.  The  assignment  of  additional  model  parameters  <P  velocity,  a  = 
6.08  km/sec;  S  velocity,  3  =  3.5  km  sec;  density,  p  -  2.74  gmcmJ:  and  Futterman 
attenuation  operator,  t =  4.5  sec  i  strongly  affects  the  moment  required  to  match 
the  observed  SH  amplitudes.  These  parameters  are  poorly  constrained  due  to 
uncertainty  regarding  the  appropriate  source  structure  and  attenuation  for  the 
mantle.  The  input  parameters  for  the  SH  synthetics  differ  from  those  used  by 
Burger  et  al.  (1986b)  in  order  to  be  consistent  with  the  EU2  structure  'Lerner- 
Lam  and  Jordan,  1987)  adopted  later  in  the  more  involved  WKBJ  spectral  SV 
modeling.  The  moments  in  Figure  2  are  the  values  required  to  match  the  observed 
peak-to-peak  27  Octobe  r  1973  SH  amplitudes  with  the  strike-slip  solution  (strike 
<b  =  287°;  dip,  5  =  90:;  rake,  K  =  0°).  Twenty-seven  nonnodal  stations  yield  an 
average  double-couple  moment  of  3.2  x  1024  dyne-cm  for  27  October  1973  and 
1.7  x  1024  dyne-cm  (20  nonnodal  observations'  for  the  2  November  1974  event. 


TABLE  1 

Event  Informa tt  n* 


Date 

UTC 

Latitude  <  .V 
Longitude  1  E 1 

777  A 

Ms 

m.  yield  'it 
Ms  yield  i 
Average  yield  ■  kt > 


27  October  197:! 
06:59:57 
70  78 
54.18 

6.94  *  0  018 
5.51  *  0  052 

4050 

28*0 

3450  r  610 


2  November  1974 
04  59  57 
70.62 
54  18 

6.78  *  0  020 
5.29  *  0  040 

2100 

1690 

1890  *  210 

Magnitudes 


*  Locations  are  from  Dahlman  and  Israelson  <  1 977 > 
and  yields  are  from  Lykes  and  Wiggins  <  1986i 
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Fig.  2.  Lone-period  SH  observations  and  half  space  synthetics  for  the  27  October  1973  explosion. 
Numbers  are  the  moments  required  to  match  tlv  observed  amplitudes  at  each  station  (  x  1024  dyne- 
cmi.  The  lower  hemisphere,  equal-area  projection  shows  the  SH  first  motions  (filled  circles  are  clock¬ 
wise*  and  the  SH  radiation  nodes  for  the  vertical  strike-slip  mechanism. 


These  moments  are  both  twice  as  large  as  the  values  found  by  Burger  et  al. 
(1986b),  due  primarily  to  the  faster  source  velocity  used  in  our  synthetics.  The 
tectonic  release  moment  for  the  1973  event  is  1.9  times  larger  than  that  for  the 
1974  explosion.  The  same  amplitude  ratio  is  also  found  for  the  long-period  P  and 
SV  amplitudes  at  common  stations  for  the  two  events  (Burger  et  al.,  1986b), 

In  order  to  decide  which  mechanism  is  correct,  it  is  necessary  to  model  the  SV 
data.  The  observed  SV  waveforms  are  remarkably  similar  for  both  events,  even 
60  sec  after  the  initial  S  arrival.  This  is  demonstrated  in  Figure  3  where  selected 
radial  and  vertical  SV  observations  from  each  event  are  compared.  The  waveform 
similarity  further  establishes  the  consistency  of  the  tectonic  release  orientation 
for  both  explosions  as  well  as  the  stability  of  the  SV  wave  Field.  Note  that  the 
SV  energy  is  distributed  over  a  long  time  window,  unlike  the  simple  SH  pulses 
in  Figure  2,  and  that  the  First  arrival  waveforms  vary  substantially.  The  azi¬ 
muthal  dependence  of  the  direct  S  arrival  is  masked  by  variation  in  the  long- 
period  SPL  amplitudes,  which  are  primarily  attributable  to  variations  in  receiver 
structure  and  epicentral  distance.  Note  also  the  prograde  elliptical  polarization 
of  the  SPL  phase  (Figure  3).  A  more  complete  azimuthal  coverage  of  the  observed 
radial  SV  data  for  the  27  October  1973  event  is  presented  in  Figure  4.  The  am¬ 
plitudes  shown  are  peak-to-peak  measurements  of  direct  S  in  millimicrons  (mpl 
normalized  to  50°  distance  by  a  geometric  spreading  correction  (Langston  and 
Helmberger,  1975).  The  direct  S  amplitudes  display  a  clear  two-lobed  azimuthal 
pattern,  with  weak  arrivals  near  152°  and  332°  (NIL,  AAM)  and  large  amplitudes 
near  62°  and  242°  (COL,  COP).  The  amplitude  variation  of  the  initial  S  arrival 
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Fig.  3.  Long-period  radial  ( R>  and  vertical  <Z)  component  SV  observations  for  27  October  1973 
and  2  November  1974  explosions.  Downward  motion  on  the  radial  components  indicates  motion  toward 
the  source.  The  numbers  are  direct  S  peak-to-peak  amplitudes  in  millimicrons  corrected  for  geometric 
spreading  to  50°  distance.  Arrows  indicate  the  arrival  time  of  the  direct  S  pulse  The  long-period 
oscillations  later  in  the  waveform  are  the  SPL  phase. 


cannot  be  attnouted  to  variations  of  receiver  structure,  and  instead,  must  rep¬ 
resent  interference  between  the  explosion  and  tectonic  release  signals.  Burger 
et  al.  (1986b,  Figure  4>  show  a  corresponding  figure  with  numerous  radial  SV 
signals  for  the  2  November  1974  explosion,  which  have  very  similar  wave  shapes. 

While  recognizing  the  limitations  of  ray  theory  for  SV  signals,  we  generated 
SV  half-space  synthetics,  including  the  explosion  i pS<  and  tectonic  release  ar¬ 
rivals  (S,  sS,  pS)  to  identify  first-order  interference  effects  These  synthetics  are 
constructed  with  the  same  parameters  used  in  the  SH  modeling.  The  explosion 
is  represented  by  a  modified  Haskell  bomb  source  (Helmberger  and  Hadley,  1981 1 
constructed  using  an  overshoot  parameter  B  =  0.1  and  rise  time  parameter  K 
—  7.5  sec  ' l.  The  long-period  explosion  SV  synthetics  are  not  especially  sensitive 
to  the  choice  of  h  and  K ;  however,  these  values  are  in  agreement  with  the  P  wave 
modeling  by  Burger  ct  al.  (1986b)  for  the  same  two  events. 

The  polarity  of  the  .  idial  component  of  direct  S  for  the  normal  fault  solution 
is  toward  the  source  at  all  azimuths,  the  same  as  for  the  explosion  (pS)  pulse. 
Accordingly,  the  two  signals  const  ructi  vely  interfere  in  the  first  cycle  of  the  direct 
S  arrival,  and  the  amplitude  of  the  smallest  observed  SV  waves  near  azimuths 
of  152°  and  332°  must  represent  the  size  of  the  explosion  signal  alone.  Thus,  the 
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Fig.  4  Long-period  radial  component  SV  observations  for  27  October  1973.  The  numbers  are  the 
observed  peak-to-peak  amplitudes  of  direct  5  in  millimicrons  normalized  to  50°  distance  using  a 
geometric  spreading  correction.  The  lower  hemisphere,  equal-area  projection  shows  the  location  of 
each  observation  and  the  SV  radiation  nodes  for  the  strike-slip  mechanism  The  large  amplitude  at 
COP  (A  =  23°)  results  from  an  upper  mantle  triplication. 


explosion  source  strength,  iK,  for  the  45°  normal  fault  solution  must  be  negligible 
to  produce  the  observed  radiation  pattern  which  has  clear  nodes.  Due  to  this, 
and  additional  problems  discussed  later  with  the  WKBJ  spectral  modeling,  the 
normal  solution  is  deemed  inappropriate. 

Ray  theory  SV  waveforms  for  the  explosion  and  the  two  remaining  plausible 
double-couple  orientations,  are  shown  separately  and  superimposed  in  Figure  5. 
In  each  case,  the  explosion  pS  amplitude  is  adjusted  to  achieve  the  desired  two- 
lobed  radiation  pattern  and  to  match  the  observed  27  October  1973  amplitudes 
at  50°  distance.  All  amplitudes  are  1.9  times  smaller  for  the  2  November  1974 
event.  Both  double-couple  orientations,  when  combined  with  the  explosion  source, 
can  produce  the  basic  amplitude  variation  observed  in  the  data;  however,  the 
azimuinal  interference  behavior  differs  for  the  two  mechanisms.  The  shape  of 
the  Green’s  function  for  the  thrust  solution  is  very  similiar  to  the  strike-slip 
solution  in  the  destructive  ( upward  direct  S )  lobe  <<b  =  242°).  The  Green's  function 
for  the  strike-slip  solution  in  the  constructive  (downward  direct  S)  lobe  (<i>  = 
332°)  is  also  similar,  but  with  opposite  polarity.  One  might  expect  that  the  simi¬ 
larity  of  the  double-couple  waveforms  over  two  out  of  four  quadrants  would  com- 
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Fic.  5.  Half-space  radial  component  S  V  synthetics  at  50’  distance  showing  first-order  interference 
effects  between  the  explosion  signal  and  the  45’-dipping  thrust  defti  and  vertical  strike-slip  i right' 
tetonic  release  solutions  The  focal  mechanism  parameters  (strike.  <t>;  rake.  h,  dip,  5i  and  azimuths 
are  appropriate  for  the  Southern  Novaya  Zemlva  test  site  The  top  row  depicts  the  explosion  waveform 
alone  with  the  peak-to-peak  amplitude  in  millimicrons  required  to  produce  the  observed  SV  radiation 
pattern  for  the  27  October  1973  event.  The  double-couple  synthetics  at  various  azimuths  are  shown 
in  the  middle  row,  and  the  waveforms  resulting  from  superposition  of  the  explosion  and  double-couple 
radiation  are  shown  in  the  bottom  row. 


plicate  efforts  to  distinguish  between  the  two  solutions,  regardless  of  the  presence 
of  additional  SV  phases.  Indeed,  for  stations  in  the  two  destructive  strike-slip 
quadrants,  no  preference  can  be  made  for  either  choice  of  focal  mechanism,  even 
for  more  complete  SV  synthetics.  However,  as  will  be  shown,  the  observations 
in  the  constructive  lobes  of  the  strike-slip  solution  allow  identification  of  the  true 
tectonic  release  orientation. 

The  thrust  double-couple  solution  (d>  =  242°,  5  =  45°,  \  =  90°)  for  the  half¬ 
space  (Mo  =  6.5  x  1024  dyne-cm,  iK  =  8.7  x  10u  cm3)  has  an  upward  first 
motion  over  all  azimuths,  opposite  that  of  the  explosion  signal,  resulting  in 
dominantly  destructive  interference.  The  slight  differences  in  the  thrust  and 
explosion  waveforms  result  in  destructive  interference  behaving  less  efficiently 
than  constructive  interference.  However,  this  interference  is  very  sensitive  to 
the  time  alignment  of  the  explosion  and  double-couple  arrivals,  which  ;s  poorly 
known.  Small  time  shifts  can  cause  much  smaller  amplitudes  at  an  azimuth  of 
242°. 

The  SV  waves  for  the  strike-slip  solution  ( M0  =  3.2  x  1024  dyne -cm,  <]/,  =  3.2 
x  10u  cm3)  have  a  sin  26  radiation  pattern.  This  results  in  alternating  lobes  of 
constructive  and  destructive  interference  in  the  first  cycle  of  the  combined  ex¬ 
plosion  plus  double-couple  synthetics.  Figure  5  shows  the  relevant  half-space 
synthetics  for  the  strike-slip  solution.  The  for  the  strike-slip  case  is  only  37 
per  cent  of  the  thrust  case,  indicating  the  importance  of  determining  the  correct 
focal  mechanism.  The  explosion  alone  must  contribute  the  energy  observed  at 
the  strike-slip  nodes,  thus,  ib*  is  constrained  to  generate  an  amplitude  of  -1900 
mti,  corresponding  to  t  e  average  observed  value  at  the  strike-slip  nodal  azi¬ 
muths.  The  synthetic  SV  amplitude  pattern  for  the  strike-slip  case  matches  the 
data  reasonably  well  for  both  events,  as  shown  in  Figure  6.  The  asymmetry  in 
the  SV  amplitude  pattern  is  caused  by  a  slight  reversal  of  the  first  motion  in  the 
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SH  SV 


Fig.  6.  Observed  and  half-space  synthetic  SH  and  SV  amplitude  patterns  for  the  27  October  1973 
and  2  November  1974  events  corrected  for  geometric  spreading  to  a  distance  of  50°.  The  source  model 
for  the  synthetic  amplitudes  is  the  explosion  plus  strike-slip  tectonic  release  model. 


destructive  quadrants.  The  synthetic  amplitudes  of  the  explosion  plus  thrust 
solution  do  not  fit  the  observed  pattern  as  well,  but  can  be  made  to  do  so  by  slight 
adjustment  ( =  1  sec)  of  the  relative  timing  of  the  arrivals.  Thus,  for  the  ray  theory 
SV  synthetics,  either  solution  is  equally  viable.  The  synthetic  SH  amplitudes 
and  polarities  over  the  full  azimuthal  range  also  closely  match  the  observations 
and  are  presented  in  Figure  6.  The  large  number  of  observations  for  both  events 
narrowly  constrain  the  possible  focal  mechanisms.  It  is  evident  from  the  simi¬ 
larity  of  interference  that  both  events  had  a  tectonic  release  component  with  the 
same  orientation. 

The  SV  ray  theory  modeling  has  eliminated  the  normal  fault  tectonic  release 
orientation;  however,  it  does  not  discriminate  between  the  thrust  and  strike-slip 
solutions.  Clearly,  the  radial  SV  waveforms  in  Figures  3  and  4  do  not  resemble 
the  simple  synthetics  of  Figure  5  for  either  solution.  It  is  necessary'  to  enlarge 
the  time  window  under  consideration  and  include  the  additional  phases  of  SV 
wave  propagation  to  reveal  any  diagnostic  interference  behavior  of  the  two  pos¬ 
sible  focal  mechanisms.  To  accomplish  this,  the  WKBJ  spectral  method  may  be 
applied  to  the  SV  propagation  problem. 

WKBJ  Spectral  Modeling 

Baag  and  Langston  (1985b)  coupled  WKBJ  theory  for  a  generalized  ray  with 
the  Watson  compound-matrix  method  to  generate  a  computationally  efficient, 
accurate  treatment  of  SV  wave  propagation.  The  method  is  based  on  Frazer 
(1977),  but  also  includes  anelasticity,  parameterized  as  complex-valued  velocity, 
and  near-source  reverberations  that  contribute  downgoing  S  energy,  particularly 
for  shallow  sources.  A  depth  is  chosen  for  coupling  the  propagator  matrix  with 
the  generalized  rays,  so  as  to  include  the  upper  mantle  regions  with  strong  ve¬ 
locity  gradients  and  low-velocity  zones  in  the  near-surface  layers.  Since  SPL 
spends  a  large  amount  of  time  in  the  crustal  wave  guide  through  multiple  re¬ 
flections  and  P-S  conversions,  the  choice  of  velocity  model  is  important.  Diffracted 
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Sp  is  likewise  sensitive  to  crustal  structure.  The  earth  model  is  discretized  above 
the  coupling  boundary  into  homogeneous  layers  and  smoothed  below  to  form  a 
continuous,  vertically  inhomogeneous  structure.  The  calculated  downgoing  re¬ 
flectivity  potential  is  used  as  the  impulse  for  the  WKBJ  treatment.  The  range 
of  ray  parameter  is  chosen  so  as  to  include  only  those  rays  that  turn  below  the 
coupled  boundary  in  the  WKBJ  domain.  The  upgoing  WKBJ  ray  is  eventually 
converted  to  the  propagator  matrix  coefficient  and  evaluated  in  the  near-receiver 
layers.  An  earth  flattening  transformation  of  the  earth  structure  is  used  for 
treatment  in  a  cylindrical  coordinate  system  (Chapman,  1973;  Muller,  1977),  and 
the  solution  is  evaluated  in  the  frequency  domain  over  a  designated  range  of 
real-valued  wavenumbers.  The  resulting  Green’s  function  is  then  convolved  with 
an  appropriate  source-time  function  and  instrument  response,  and  inverse  Four¬ 
ier  transformed  to  produce  displacements  in  the  time  domain. 

The  WKBJ  spectral  method  requires  specification  of  numerous  parameters  to 
constrain  the  earth  model,  source  iype,  and  various  computational  options.  Three 
earth  models  were  considered  in  this  study:  the  EU2  model  of  Lerner-Lam  and 
Jordan  (1987),  and  models  TNA  and  SNA  of  Grand  and  Helmberger  (1984).  These 
three  structures  are  shown  to  600  km  depth  in  Figure  7.  Below  this  depth,  the 
three  models  have  the  same  shear  velocity.  The  TNA  and  SNA  structures  are 
extended  down  to  the  core-mantle  boundary  with  the  Jeffreys-Bullen-A  (JBA) 
model.  Since  it  is  not  necessary  to  compute  P  potentials  deep  in  the  mantle  for 
teleseismic  SV  waves,  the  P-wave  velocity  is  only  defined  in  the  homogeneous 
near-surface  layers.  For  the  TNA  and  SNA  models,  the  P-wave  velocity  structure 
is  prescribed  by  the  JBA  model.  The  SV  wave  train  is  sensitive  to  crustal  structure 
and  each  of  the  three  earth  models  promote  unique  behavior;  thus,  it  was  desired 
to  use  one  model  for  the  majority  of  the  modeling  in  order  to  isolate  structural 
effects  from  other  variables.  As  shown  in  Figure  7,  the  EU2  model  generally  has 
velocities  intermediate  between  TNA  and  SNA.  Additionally,  EU2  was  formu¬ 
lated  by  a  waveform  inversion  of  fundamental  and  higher  mode  Rayleigh  waves 
propagating  across  northern  Eurasia.  Accordingly,  it  seems  appropriate  to  use 
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EU2  for  the  Novaya  Zemlya  source  structure  throughout  the  modeling,  and  when¬ 
ever  possible  for  the  receiver  structure. 

A  Q  structure  similar  to  SL8  (Anderson  and  Hart,  1978)  was  used  in  all  syn¬ 
thetics  (Figure  7  >.  The  effects  of  the  frequency-dependent  Q  model  (EURS)  of  Der 
et  al.  (1986),  formulated  from  seismic  observations  for  Eurasian  seismic  paths, 
are  addressed  later  in  the  "Discussion.”  The  boundary  between  the  homogeneous 
layers  and  the  vertically  inhomogeneous  half-space  is  set  at  200  km  depth  to 
keep  the  WKBJ  ray  turning  smoothly  below  the  strong  velocity  gradients  and 
low-velocity  zone  in  the  crust.  The  source  depth  is  placed  at  1  km  for  both  the 
explosion  and  the  tectonic  release  component.  Although  the  depth  of  the  double¬ 
couple  source  is  not  constrained  as  well  as  for  the  explosion,  there  is  no  evidence 
for  distinct  depth  phase  arrivals  in  either  the  SH  or  SV  data.  Depths  from  0.5 
to  1.5  km  were  also  explored,  but  the  effect  on  the  waveforms  and  amplitudes  is 
negligible.  This  is  not  particularly  surprising  considering  the  limited  resolution 
of  the  synthetic  Green’s  functions,  which  are  calculated  with  128  frequency  points 
and  a  2.0-sec  sampling  interval.  This  coarse  sampling  rate  was  found  to  be  ad¬ 
equate  for  our  comparisons  with  long-period  data.  The  higher  frequency  infor¬ 
mation  is  not  crucial  to  our  understanding  of  the  general  wave  behavior,  and  the 
additional  computation  time  required  for  a  smaller  sampling  interval  is  not  jus¬ 
tified.  Because  the  convolution  of  the  source-time  function  and  instrument  re¬ 
sponse  is  discretized  at  a  2.0-sec  interval,  longer  duration  source-time  functions 
are  preferred  to  avoid  sampling  errors.  As  before,  a  3.0-sec  trapezoid  is  used  for 
the  tectonic  release  component  source-time  function,  and  a  modified  Haskell 
bomb  source  with  B  =  0.1,  and  K  =  7.5  sec1  describes  the  isotropic  source 
component.  These  source-time  functions  are  consistent  with  results  from  relative 
waveform  analysis  by  Burger  et  al.  (1986a)  for  the  Southern  Novaya  Zemlya 
explosions.  Although  the  sensitivity  of  the  synthetics  to  the  different  source-time 
functions  is  limited,  care  has  been  exercised  to  maintain  the  absolute  amplitude 
information  for  each  source. 

For  the  purposes  of  this  study,  we  have  restricted  our  attention  to  observations 
in  the  distance  range  of  30°  to  60°.  The  different  sources  considered  here  are  most 
distinct  in  this  distance  range.  Beyond  60°,  the  SV  phases  other  than  direct  S 
become  less  important  and  the  SV  waveforms  can  be  accurately  represented  with 
ray  theory  modeling  (Langston  and  Baag,  1985).  Furthermore,  receivers  between 
30°  and  60°  are  not  influenced  by  upper  mantle  triplications  or  ScS  and  SKS 
phases.  Vertical  and  radial  components  of  the  WKBJ  spectral  synthetics  for  both 
the  explosion  and  double-couple  sources  are  displayed  in  Figure  8  for  the  EU2 
source  and  receiver  structure.  In  Figure  8,  the  double-couple  waveforms  shown 
correspond  to  either  the  thrust  signals  at  all  azimuths  or  the  destructive  quad¬ 
rants  of  the  strike-slip  solution.  Recall  that  the  shape  of  the  Green's  function  for 
the  thrust  solution  at  all  azimuths  is  very  similar  to  that  of  the  strike-slip  solution 
in  the  two  destructive  quadrants.  The  Green’s  function  for  the  strike-slip  solution 
in  the  constructive  lobe  also  has  the  same  shape,  but  with  opposite  polarity. 

The  important  phases  of  the  SV  wave  train  are  identified  in  Figure  8.  The  Sp 
conversion  immediately  precedes  S  as  a  small  amplitude  arrival  of  opposite  po¬ 
larity  on  the  radial  component.  The  exact  timing  and  strength  of  the  Sp  arrival 
are  dependent  on  crustal  velocities  and  depths  to  velocity  discontinuities  under 
the  receiver.  This  phenomenon  of  Sp  precursors  to  direct  S  has  been  investigated 
by  Jordan  and  Frazer  (1975)  and  Burdick  and  Langston  (1977)  as  a  means  of 
discerning  receiver  structure.  The  long-period  SPL  phase  is  frequently  the  most 
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Fig.  8.  Vertical  and  radial  SV  WKBJ  spectral  synthetics  for  the  explosion  and  double-couple 
sources  in  the  distance  range  30“  to  60'.  The  aouble-couple  waveforms  correspond  to  the  thrust  signals 
at  all  azimuths  and  the  destructive  quadrants  of  the  strike-slip  solution  The  important  phases  of  the 
SV  wave  train  are  identified  except  for  diffracted  Sp.  which  arrives  after  direct  S  at  distances  greater 
than  50°. 


diagnostic  feature  of  the  SV  wave  train.  The  normal  dispersion  of  the  phase  is 
evident  over  the  30°  range  shown,  and  the  SPL  amplitudes  gradually  decrease 
with  greater  distance.  SPL  is  preferentially  excited  by  the  double-couple  source 
relative  to  the  explosion  source,  particularly  at  closer  distances.  Additionally, 
the  overall  synthetic  explosion  waveform  contains  more  high  frequency  energy 
than  the  double-couple  source.  Even  if  the  Sp  precursor  is  ignored,  the  synthetics 
in  Figure  6  are  considerably  different  in  the  first  20  sec  from  the  half-space 
synthetics  shown  in  Figure  5.  The  rigorous  treatment  of  the  SV  propagation 
problem  afforded  by  the  WKBJ  spectral  method  produces  synthetic  waveforms 
that  begin  to  include  the  complexities  seen  in  Figures  3  and  4.  The  accurate 
computation  of  the  Green’s  function  becomes  even  more  crucial  when  the  two 
sources  are  superimposed.  No  time  shift  of  the  double-couple  relative  to  the  ex¬ 
plosion  source  is  utilized  for  any  of  the  modeling,  but  small  time  shifts  prove  less 
important  than  in  the  half-space  modeling  because  we  are  comparing  the  full 
waveform. 

The  two  focal  mechanisms  combined  with  an  appropriate  size  explosion  source 
are  shown  in  Figure  9  over  a  90°  azimuthal  range.  In  comparison  to  the  ray  theory 
amplitudes,  the  WKBJ  spectral  method  amplitudes  for  the  thrust  solution  (M0 
=  6.5  x  1024  dyne-cm,  di*  =  3.8  r  0.6  x  1011  cm3)  better  match  the  variation 
seen  in  the  data  (Figure  6),  and  the  d),  required  to  accomplish  this  is  less  than 
half  of  the  half-space  value.  This  change  in  di«  can  be  attributed  to  components 
in  the  WKBJ  spectral  waveforms  producing  slightly  more  efficient  destructive 
interference  effects  than  are  possible  with  the  half-space  modeling.  The  strike- 
slip  solution  ( Af0  =  3.2  x  1024  dyne-cm,  di*  =  3.8  *  0.5  x  10“  cm3)  also  matches 
the  amplitude  variation  of  Figure  6,  but  di*  is  roughly  the  same  as  for  the  half¬ 
space  modeling. 

The  synthetic  seismograms  of  Figures  8  and  9  were  constructed  using  an  EU2 
earth  structure.  The  suitability  of  EU2  for  the  source  structure  has  been  dis¬ 
cussed.  As  a  receiver  structure,  EU2  performs  satisfactorily  in  most  cases,  and 
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a  constant  receiver  structure  is  preferred  to  isolate  the  interference  effects  of 
each  solution.  However,  the  relationship  between  receiver  structure  and  the  re¬ 
sulting  SV  phases  is  quite  complicated,  and  it  is  therefore  necessary  to  examine 
different  source-receiver  paths  to  establish  the  structural  effects.  Figure  10  de¬ 
picts  WKBJ  spectral  seismograms  for  the  strike-slip  solution  combined  with  the 
explosion  at  40°  distance  for  numerous  source-receiver  structural  combinations. 
Throughout  this  modeling,  the  source  depth  is  maintained  at  1  km,  and  as  a 


Fig.  9.  Radial  SV  WKBJ  spectral  synthetics  for  each  tectonic  release  solution  combined  with  an 
appropriate  size  explosion  over  a  representative  90°  azimuthal  range.  The  45°-dipping  thrust  signal 
interferes  destructively  with  the  explosion  at  all  azimuths.  The  four-lobed  strike-slip  solution  signal 
alternates  between  destructive  and  constructive  interference  in  opposite  lobes.  The  amplitudes  are 
in  millimicrons  and  are  appropriate  for  the  27  October  1973  event. 
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Fig.  10.  Radial  SV  WKBJ  spectral  synthetics  for  the  strike-slip  solution  superimposed  with  the 
explosion  signal  for  various  source-receiver  structural  combinations  (E.  EU2;  T,  TNA;  S,  SNA).  The 
waveforms  are  equivalent  when  source  and  receiver  structures  are  exchanged  due  to  near-reciprocity 
with  a  1  km  source  depth.  Numbers  are  peak-to-peak  amplitudes  in  millimicrons  for  a  distance  of 
40°. 
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result,  near-reciprocity  exists  when  source  and  receiver  structures  are  inter¬ 
changed.  Furthermore,  this  near-reciprocity  indicates  that  source  and  receiver 
structures  are  equally  important  in  the  resultant  seismogram  for  shallow'  sources. 
The  waveforms  in  Figure  10  differ  most  dramatically  in  the  relative  excitation 
of  the  SPL  phase,  the  timing  of  the  initial  arrival,  and  the  peak-to-peak  direct 
S  amplitudes.  The  variance  in  amplitude  and  arrival  time  is  an  expected  result 
of  the  different  velocities  of  each  structural  combination.  The  amplitude  infor¬ 
mation  in  Figure  10  does  not  help  constrain  the  appropriate  earth  structure,  since 
source  strengths  may  be  weighted  accordingly  to  fit  the  observations.  Likewise, 
the  timing  of  the  initial  S  arrival  is  not  resolved  well  enough  in  the  observations 
for  diagnostic  analysis.  Relative  excitation  of  the  SPL  phase,  however,  may  be 
used  to  distinguish  earth  structure.  With  an  EU2  source  structure,  the  EU2  and 
SNA  receiver  structures  yield  nearly  identical  wave  shapes  and  amplitudes.  The 
TNA  receiver  structure  increases  SPL  relative  to  direct  S  more  than  SNA  or 
EU2  for  the  distance  range  considered.  This  result  differs  from  that  of  Baag  and 
Langston  (1985a)  who  found  SNA  excited  SPL  more  than  TNA.  This  apparent 
inconsistency  results  from  the  longer  period  30-100  (To  =  30  sec.  tg  =  100 


sec)  WWSSN  instrument  and  long-duration 


1 

t 


source-time 


function  used  in  their  synthetics. 

An  additional  complication  for  modeling  the  Novaya  Zemlya  data  is  the  as¬ 
sumption  of  lateral  homogeneity  in  earth  structure.  Evidence  has  been  presented 
for  the  existence  of  strong  lateral  variations  in  upper  mantle  velocities  underlying 
western  Eurasia  (England  et  al.,  1977;  Rial  et  al.,  1984).  Lateral  heterogeneity 
within  the  source  or  receiver  structure  could  perturb  the  SV  waveform  in  a  man¬ 
ner  that  is  difficult  to  identify.  Our  modeling  allows  some  flexibility  in  receiver 
structure  that  may  in  some  cases  account  for  lateral  heterogeneity  and  other 
structural  complexities.  With  the  availability  of  representative  Green's  functions 
for  the  explosion  and  double-couple  sources,  and  insight  into  the  appropriate 
earth  model  and  relative  weighting  necessary  for  each  source  to  produce  the 
observed  amplitude  variations,  the  task  of  isolating  a  preferred  double-couple 
orientation  may  be  approached. 


Comparison  of  Synthetic  and  Observed  Waveforms 

The  strike-slip  and  thrust  tectonic  release  orientations  are  combined  with  an 
appropriate  size  explosion  source  to  produce  the  synthetic  waveforms  compared 
with  the  27  October  1973  observations  in  Figure  11.  The  uniform  shape  of  the 
two  synthetic  seismograms  at  POO  and  BLC  are  typical  of  a  number  of  stations 
in  the  strike-slip  destructive  quadrants.  For  the  90°  azimuth  range  spanned  by 
each  of  these  lobes,  the  interference  of  both  tectonic  release  solutions  is  destruc¬ 
tive,  and  the  resulting  waveforms  are  typically  quite  similar,  suppressing  any 
preference  for  either  solution.  However,  the  observed  waveforms  at  all  of  the 
other  stations  are  c' early  matched  better  by  the  strike-slip  solution.  In  fact,  the 
azimuthal  range  where  the  strike-slip  solution  is  preferred  actually  extends  out 
of  the  180°  range  that  defines  the  constructive  lobes.  Stations  in  the  destructive 
quadrants  located  near  tno  strike-slip  nodes  (INK)  contain  relative  source  con¬ 
tributions  that  are  still  considerably  different  from  the  thrust  solution.  For  all 
stations  where  the  two  double-couple  solutions  produce  distinct  waveforms,  the 
strike-slip  orientation  consistently  produces  superior  matches  with  the  observed 
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I  60  sec  1 


Fig.  11.  Radial  SV  observations  for  27  October  1973  compared  with  the  WKBJ  spectral  synthetic 
waveforms  for  the  45°  thrust  and  vertical  strike-slip  solutions  combined  with  the  explosion  source. 
P- wave  radiation  nodes  for  the  two  solutions  are  shown  in  the  lower  hemisphere  projection.  The  source- 
receiver  structures  follow  the  notation  of  Figure  12.  Numbers  shown  are  the  observed  and  synthetic 
peak-to-peak  amplitudes  in  millimicrons. 


seismograms.  Given  the  complexities  intrinsic  to  the  SV  waveform  for  both  dou¬ 
ble-couple  sources  and  the  explosion,  it  is  remarkable  how  well  the  synthetic 
seismograms  resemble  the  data.  The  quality  of  the  matches,  and  the  uniqueness 
and  preference  of  the  strike-slip  solution  over  a  considerable  azimuthal  range, 
permit  the  elimination  of  the  thrust  orientation  of  tectonic  release  for  the  two 
explosions  considered.  Thus,  careful  forward  modeling  of  SV  wave  propagation 
can  resolve  the  appropriate  tectonic  release  double-couple  in  situations  where 
other  phases  and  ray  theory  modeling  fail  to  yield  a  unique  solution. 

A  large  suite  of  synthetic  SV  seismograms  are  compared  with  observed  wave¬ 
forms  in  Figures  12  and  13,  for  the  27  October  1973  and  2  November  1974  events, 
respectively.  The  synthetics  are  constructed  using  the  aforementioned  input  pa¬ 
rameters  and  the  strike-slip  solution.  The  receiver  structure  is  indicated  by  a 
star  if  different  from  EU2.  The  traces  have  not  been  filtered;  hence,  the  observed 
(top)  seismograms  contain  slightly  more  high-frequency  energy  than  the  2.0-sec 
sampled  synthetic  (bottom)  traces.  Nonetheless,  the  two  signals  match  remark¬ 
ably  well. 

Gradations  in  relative  source  contribution  and  changing  polarity  cause  the 
waveforms  to  vary  considerably  over  azimuth.  Most  of  the  existing  misfit  is  read¬ 
ily  attributable  to  structural  effects;  misrepresentation  of  the  actual  receiver 
structure  by  the  idealized  models.  In  some  cases,  the  first-order  differences  can 
be  diminished  by  appealing  to  either  TNA  or  SNA  for  the  receiver  structure. 
This  is  only  necessary  for  the  few  anomalous  stations  indicated.  The  effects  of 
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Fig  12.  Observed  SV  seismograms  for  27  October  1973  compared  with  WKB.I  spectral  synthetics 
using  the  strike-slip  tectonic  release  solution.  Stations  are  displayed  in  order  of  increasing  azimuth 
Numbers  shown  are  the  observed  and  synthetic  peak-to-peak  amplitudes  in  millimicrons  Source- 
receiver  structures  are  indicated  in  the  key 


each  receiver  structure  for  the  EU2  source  model  has  been  previously  addressed 
and  shown  in  Figure  10.  We  have  not  attempted  the  difficult  task  of  perturbing 
the  earth  models  to  improve  the  matches.  The  good  quality  of  the  synthetic- 
observed  fits  suggests  that  if  receiver  structures  are  well-known  for  a  number 
of  stations,  it  may  be  possible  to  modify  the  EU2  source  model  to  better  under¬ 
stand  the  seismic  structure  underlying  Novaya  Zemlya.  An  additional  source  of 
waveform  mismatch  arises  from  the  modeling  assumption  of  horizontal  velocity 
discontinuities.  Dipping  layers,  notably  a  dipping  Moho.  may  alter  the  dominant 
period  and  amplitude  of  the  Sp  and  SPL  phases.  The  behavior  of  SPL  at  MAT, 
1ST,  and  TOL  is  likely  due  to  receiver  complexities  not  accounted  for  in  the  plane- 
layered  earth  models.  Hence,  the  choice  of  receiver  model  is  not  intended  to 
directly  reveal  the  actual  receiver  velocity  structure.  Rather,  it  provides  a  mod- 
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Fig.  13.  Same  as  Figure  12.  but  for  the  2  November  1974  event. 


eling  flexibility  that  accounts  to  first-order  for  lateral  heterogeneity,  variable 
low-velocity  zones,  and  other  indeterminate  structural  complexities. 

The  peak-to-peak  amplitudes  in  mp  for  the  initial  S  arrivals  are  shown  in 
Figures  12  and  13.  For  most  stations,  the  amplitudes  agree  within  a  factor  of  2. 
This  is  to  be  expected,  given  that  the  observed  amplitude  pattern  was  initially 
used  to  constrain  the  relative  source  contributions.  The  success  of  this  application 
of  the  WKBJ  spectral  method  encourages  application  of  this  procedure  to  related 
problems,  chiefly  because  it  has  involved  the  calculation,  relative  weighting,  and 
combination  of  two  very  complex  Green’s  functions. 

F  Factors 

The  F  factor  was  introduced  by  Toksdz  et  al.  ( 1965)  as  a  means  of  characterizing 
the  ratio  of  nonisotropic  to  isotropic  energy, 

F  =  \.5(ModciMo'’xp) 

where  the  moment  of  the  double-couple  may  be  determined  by  SH  modeling  and 
the  explosion  moment  by. 


Moexp  —  4Ttpa2uj» 

(Aki  et  al.,  1974).  The  determination  of  Mvxp  depends  primarily  on  the  seismic 
velocity  and  attenuation  structure  assumed  in  the  modeling.  The  poor  resolution 
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of  these  parameters  introduces  a  baseline  uncertainty  lor  the  absolute  Mo  values. 
Although  the  moments  are  not  well  constrained,  the  ratio  of  energy  contributions 
should  remain  nearly  constant,  regardless  of  the  earth  model  used. 

For  the  assumed  source  values  of  a  (6.08  km/sec)  and  p  (2.74  gmcm3),  the 
double-couple  moment  is  3.2  x  1024  dyne-cm  for  the  27  October  1973  event  and 
1.7  x  1024  dyne-cm  for  the  2  November  1974  event.  The  iK  determined  from 
strike-slip  SV  half-space  modeling  (3.2  x  1011  cm3)  translates  to  M0exs  =  4.2  x 
1024  dyne-cm,  and  F  =  1.2  for  27  October  1973.  The  M{fxp  and  M0dc  for  2  No¬ 
vember  1974  are  both  1.9  times  smaller  than  the  27  October  1973  values,  pro¬ 
ducing  the  same  F  factor.  Using  the  more  complete  WKBJ  spectral  results  (di* 
=  3.8  x  10u  cm3  for  27  October  1973  and  >1/,  =  2.0  x  io11  cm3  for  2  November 
1974),  F  =  1.0  for  27  October  1973  and  2  November  1974.  These  F  factors  are 
within  a  factor  of  2  of  the  results  found  by  Burger  et  al.  1 1986b)  using  a  combi¬ 
nation  of  long-period  P- wave  and  S- wave  modeling.  Given  the  inherent  uncer¬ 
tainties  that  limit  resolution  of  the  F  factor,  it  is  not  possible  to  argue  a  preference 
for  either  set  of  results.  However,  in  both  studies,  the  F  factor  is  relatively  small. 
Consequently,  yields  based  on  azimuthally  distributed  Ms  measurements  for 
these  events  (Sykes  and  Wiggins,  1986)  should  not  be  significantly  biased  by 
tectonic  release  contamination. 


Discussion 

An  unconstrained  deviatoric  component  accompanying  explosion  signals  can 
potentially  alter  yield  estimates  derived  from  surface  waves,  pS  amplitudes,  or 
P  waves.  Rayleigh  and  Love  wave  excitation  by  explosive  sources  has  been  in¬ 
vestigated  by  numerous  authors,  but  it  is  not  possible  to  uniquely  determine  the 
tectonic  release  mechanism  from  fundamental  modes  alone  ( Aki  and  Tsai,  1972: 
Toksoz  and  Kehrer,  1972;  Given  and  Mellman.  1986).  Patton  1 1988)  included 
higher  mode  surface  waves  to  resolve  the  deviatoric  moment  tensor,  but  this 
method  requires  knowledge  of  the  source  spall  characteristics.  SV  amplitudes 
have  been  successfully  used  in  this  study  and  others  (Lay  et  al,  1984b)  to  de¬ 
termine  explosion  source  strength  when  tectonic  release  is  either  of  known  ori¬ 
entation  or  very  small.  Helmberger  1 1985 >  used  pS  amplitudes  at  nodal  tectonic 
release  stations  to  obtain  accurate  yield  estimates  for  NTS  events.  Wallace  et  al. 
(1983)  identified  a  strong  sP  signature  in  long-period  P  waves  that  indicates  a 
high-frequency  component  of  tectonic  release  radiation.  Additionally.  Lay  et  al. 
(1984a)  examined  short-period  P  waves  and  concluded  that  it  is  possible  that 
tectonic  release  can  affect  teleseismic  short-period  amplitudes.  Assessment  of 
explosion  source  strength  using  any  of  these  phases  requires  the  determination 
of  tectonic  release  moment  and  focal  mechanism. 

Application  of  the  WKBJ  spectral  method  has  been  shown  to  be  an  effective 
tool  for  discerning  the  appropriate  tectonic  release  mechanism.  Just  as  important, 
however,  are  the  implications  for  further  discrimination  of  source  velocity  struc¬ 
ture  and  attenuation  along  the  source-receiver  path.  As  previously  mentioned, 
the  EU2  source  model  could  possibly  be  perturbed  to  produce  uniform  improve¬ 
ments  of  the  observed-synthetic  fits.  This  endeavor  would  require  detailed  in¬ 
formation  on  receiver  structures  for  a  suite  of  stations  with  little  noise  contam¬ 
ination,  which  may  be  po  sible  using  the  data  set  from  these  two  large  Novaya 
Zemlya  events.  In  addition  to  further  understanding  of  source  velocity  structure, 
an  appropriate  attenuation  model  could  perhaps  be  better  resolved  with  system¬ 
atic  WKBJ  spectral  modeling.  A  frequency-dependent  Q  model  (EURS)  was  re- 
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Fig.  14.  Radial  SV  WKBJ  spectral  synthetics  comparing  EURS  with  the  modified  SL8  attenuation 
structure  for  the  explosion  ana  strike-slip  sources  in  an  EU2  earth.  Numbers  are  peak-to-peak  am¬ 
plitudes  in  millimicrons  at  50°  distance  and  strike.  *  =  332°. 


cently  proposed  by  Der  et  al.  (19861  using  a  large  data  set  of  teleseismic  body 
waves  for  Eurasian  shield  areas.  As  shown  in  Figure  7,  the  EURS  model  has 
significantly  higher  Q  values  in  comparison  to  the  modified  SL8  Q  structure  used 
previously.  As  a  result,  using  model  EURS  produces  synthetic  waveforms  con¬ 
taining  more  high-frequency  information  and  higher  amplitudes  for  the  same 
relative  source  strengths.  Synthetic  seismograms  for  explosion  and  vertical 
strike-slip  sources  are  presented  in  Figure  14  for  an  EU2  earth  with  the  EURS 
Q  structure  at  0.03  Hz.  Although  absolute  amplitudes  will  still  be  ambiguous, 
the  frequency  content  of  synthetic  waveforms  for  a  given  Q  structure  could  be 
used  to  identify  a  Q  model  for  Novaya  Zemlya  that  is  more  appropriate  than  the 
modified  SL8  version  we  used.  Furthermore,  exploration  of  these  additional  pa¬ 
rameters  would  likely  yield  improved  matches  with  some  of  the  observed  signals, 
particularly  if  the  Green’s  functions  are  sampled  at  shorter  time  intervals  than 
used  in  this  study. 


Conclusions 

We  have  demonstrated  that  quantitative  modeling  of  complete  teleseismic  SV 
waveforms  for  explosions  with  tectonic  release  is  viable  and  can  provide  greater 
resolution  of  the  double-couple  orientation  than  modeling  of  SH  data  alone.  As¬ 
suming  that  the  nonisotropic  source  is  coincident  in  space  and  time  with  the 
explosion  source,  WKBJ  spectral  synthetic  SV  seismograms  have  been  con¬ 
structed,  including  Sp ,  diffracted  Sp ,  and  SPL  contributions.  We  apply  this  pro¬ 
cedure  to  SV  signals  from  large  underground  explosions  for  which  SH  modeling 
indicates  that  the  tectonic  release  orientation  is  either  vertical  strike-slip  or  453- 
dipping  thrust.  Comparison  of  synthetic  and  observed  SV  waveforms  for  two 
Southern  Novaya  Zemlya  events  strongly  suggests  that  a  vertical  strike-slip 
mechanism  is  appropriate  for  the  27  October  1973  and  2  November  1974  Southern 
Novaya  Zemlya  events.  In  general,  the  resolution  of  true  tectonic  release  ori¬ 
entation  has  important  ramifications  for  estimates  of  explosion  yield  using  sur¬ 
face  waves  (Af.sl,  P  waves  (mil,  or  pS  amplitudes.  Thus,  this  modeling  approach 
may  be  useful  for  other  events  for  which  we  lack  a  well-constrained  tectonic 
release  mechanism.  The  success  of  this  first-order  modeling  effort  also  suggests 
that  additional  understanding  of  source  velocity  and  attenuation  structure  is 
possible  with  systematic  WKBJ  spectral  modeling. 
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Observations  of  Teleseismic  P  Wave  Coda  for  Underground  Explosions 
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Abstract — The  early  P  wave  coda  (5-15  sec  after  the  first  arrival)  of  underground  explosions  at  the 
Nevada  Test  Site  is  studied  in  the  time  domain  using  2082  teleseismic  short-period  recordings,  with  the 
intent  of  identifying  near-source  contributions  to  the  signals  in  the  frequency  range  0.2-2.0  Hz.  Smaller 
magnitude  events  tend  to  have  relatively  high  coda  levels  in  the  0  4-0.8  Hz  frequency  band  for  both  Yucca 
Flat  and  Pahute  Mesa  explosions.  Coda  complexity  in  this  low-frequency  passband  is  negatively  correlated 
with  burial  depth  for  Pahute  Mesa  events  but  is  only  weakly  correlated  with  depth  for  Yucca  Flat  events. 
Enhanced  excitation  of  relatively  long-period  scattered  waves  for  smaller,  less  deeply  buried  events  is 
required  to  explain  this  behavior.  Coda  complexity  in  the  0  8  11  Hz  band  is  positively  correlated  with 
magnitude  and  depth  for  Pahute  Mesa  events,  but  has  no  such  dependence  for  Yucca  Flat  events.  This  may 
result  from  systematic  variations  between  the  spectra  of  direct  signals  and  coda  arrivals  caused  by  pP 
interference  for  the  largest  events,  all  of  which  were  detonated  at  Pahute  Mesa  Another  possible  explana¬ 
tion  is  a  frequency-dependent  propagation  effect  on  the  direct  signals  of  the  larger  events,  most  of  which 
were  located  in  the  center  of  the  mesa  overlying  strong  lateral  velocity  gradients  in  the  crust  and  upper 
mantle.  Event  average  complexity  varies  spatially  for  both  test  sites,  particularly  in  the  0.8-1 . 1  Hz  band, 
providing  evidence  for  frequency-dependent  focussing  or  scattering  by  near-source  structure.  Both  the 
direct  arrivals  and  the  early  coda  have  strong  azimuthal  amplitude  patterns  that  are  produced  by  defo- 
cussing  by  mantle  heterogeneity  The  direct  arrivals  have  stronger  coherent  azimuthal  patterns  than  the 
early  coda  for  Pahute  Mesa  events,  indicating  more  pronounced  deep  crustal  and  shallow  mantle  defo- 
cussing  for  the  direct  signals  However,  for  Yucca  Flat  events  the  direct  arrivals  have  less  coherent 
azimuthal  patterns  than  the  coda,  suggesting  that  a  highly  variable  component  of  near-source  scattering 
preferentially  affecting  the  downgoing  energy  is  superimposed  on  a  pattern  produced  by  mantle  hetero¬ 
geneity  that  affects  the  entire  signal.  This  complicated  behavior  of  the  direct  arrivals  may  be  the  result  of 
triplications  and  caustics  produced  by  the  complex  basement  structure  known  to  underlie  the  Yucca  Flat 
test  site.  The  presence  of  strong  azimuthal  patterns  in  the  early  coda  indicates  that  source  strength 
estimates  based  on  early  coda  are  subject  to  biases  similar  to  those  affecting  estimates  based  on  direct 
arrivals 


Key  words:  Scattering,  coda,  underground  explosions,  yield  estimation,  wave  propagation 


Introduction 

Reliable  determination  of  the  energy  released  by  seismic  sources  is  a  fundamental 
problem  in  earthquake  and  explosion  seismology.  For  relatively  small  energy  sources, 
it  is  necessary  to  utilize  short-period  seismic  waves  to  determine  the  source  strengths, 
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although  it  is  well-known  that  the  scale  lengths  of  heterogeneity  within  the  earth  are 
such  that  propagational  effects  are  particularly  pronounced  for  the  short-period 
waves  of  interest.  Thus,  there  has  been  a  concentrated  effort  to  achieve  a  better 
theoretical  and  empirical  understanding  of  short-period  wave  propagation  processes 
in  the  earth.  The  primary  objective  in  the  field  of  explosion  seismology  is  to  develop 
a  procedure  for  seismic  yield  estimation  using  short-period  waves  that  is  both  precise 
(i.e.,  consistent)  and  accurate. 

The  currently  preferred  methods  of  yield  estimation  employ  seismic  magnitude 
measurements,  which  can  be  divided  into  two  classes:  those  that  rely  on  a  direct 
phase,  such  as  nttf",  (based  on  the  first  cycle  of  the  short-period  P  waveform);  and 
those  that  are  based  on  a  scattered  phase,  such  as  Lg  or  P  coda.  The  advantage  of 
direct  phases  is  that  they  can  be  treated  to  a  large  degree  as  deterministic  signals, 
which  are  generally  expected  to  be  rather  simple  for  underground  explosions.  This 
allows  quantitative  application  of  corrections  for  predictable  biases  arising  from 
source  function  scaling,  free-surface  interaction,  attenuation  vanations  and  propaga¬ 
tion  effects  due  to  known  large-scale  three-dimensional  structure  (Lay  et  a!..  1984; 
Lay.  1985;  Cormier.  1987).  However,  the  earth  is  heterogeneous  to  such  a  degree 
that  there  appears  to  be  a  limit  to  the  precision  that  can  be  obtained  in  yield  estima¬ 
tion  using  direct  phases  (Bache.  1982).  On  the  other  hand,  scattered  phases  must 
usually  be  treated  stochastically,  but  they  do  have  the  attribute  that  the  scattered 
arrivals  sample  many  raypaths.  presumably  averaging  out  to  a  large  extent  the  vari¬ 
ous  heterogeneities  sampled  (Aki.  1982;  Cormier.  1987).  Several  authors  have  ex¬ 
amined  the  relative  performance  of  magnitudes  based  on  teleseismic  P  coda  and 
direct  P  as  relative  size  estimators,  concluding  that  P  coda  generally  yields  more 
precise  size  estimates  (  Bullitt  and  Cormier.  1984;  Baumgardt.  1985;  Gupta  et 
al. ,  1985;  Lay  and  Welc.  1987) 

In  addition,  approximately  half  of  the  energy  in  teleseismic  P  coda  from  shallow 
sources  appears  to  result  from  near-source  scattering  (  Dainty.  1985).  Thus,  the  P 
coda  may  provide  information  about  the  source  region,  either  through  use  as  a 
reference  for  examining  direct  signal  behavior  (Gupta  and  Blandford.  1987).  or 
by  analysis  of  the  coda  itself  (Lay  and  Well.  1987;  Lay.  1 987a, b).  While  some 
empirical  progress  has  been  made  in  using  P  coda  for  both  size  estimates  and  deter¬ 
mining  source  properties,  the  exact  origins  of  P  coda  remain  somewhat  enigmatic.  It 
presumably  is  generated  by  scattering  from  small-scale  heterogeneities  and  may 
therefore  be  useful  for  characterizing  crust  and  mantle  velocity  variations.  However, 
the  relative  importance  of  processes  such  as  mode  conversions,  forward  scattering  or 
back  scattering  is  not  well  established.  Improving  our  knowledge  of  the  processes 
involved  in  the  generation  of  P  coda  is  essential  to  a  basic  understanding  of  short- 
period  wave  propagation  in  the  earth. 

In  this  paper,  we  will  examine  the  properties  of  the  early  P  coda  (defined  here  as 
the  first  5-15  sec  after  the  first  arrival  of  the  P  wave)  from  underground  explosions 
at  the  Nevada  Test  Site.  The  early  coda  is  particularly  interesting  since  there  is  some 
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possibility  of  treating  it  deterministically  as  well  as  statistically,  because  the  near- 
source  contributions  must  be  generated  within  a  restricted  volume  around  the  source. 
Deterministic  treatments  that  attempt  to  explain  particular  arrivals  in  the  coda  are 
presented  by  Boucher  ( 1973),  who  analyzed  the  first  3  sec  of  near-source  record¬ 
ings,  and  by  Lay  ( 1987b)  who  inverted  for  isotropic  scatterers  using  the  first  15  sec 
of  the  P  wavetrain.  However,  inversions  of  the  latter  type  are  still  very  difficult  and 
poorly  posed  due  to  our  lack  of  knowledge  of  the  predominant  scattering  mecha¬ 
nisms  producing  the  coda;  accordingly,  in  this  study,  we  will  examine  the  statistical 
properties  of  the  early  P  coda  empirically,  to  motivate  and  help  design  future  inver¬ 
sions  for  scatterers.  These  properties  will  also  provide  an  observational  basis  for 
evaluating  the  size  estimation  capabilities  of  the  early  coda. 

The  statistical  properties  of  early  P  coda  have  already  been  extensively  studied  for 
explosions  at  the  Pahute  Mesa  test  site,  as  well  as  for  the  two  Novaya  Zemlya  test 
sites  and  the  Amchitka  test  site  by  Lay  and  Welc  ( 1987)  and  Lay  ( 1987a).  One  of 
the  interesting  results  obtained  in  these  data-intensive  studies  was  the  detection  of  a 
negative  correlation  of  coda  "complexity"  (power  in  the  coda  relative  to  the  direct 
arrival)  with  both  magnitude  and  depth  of  burial  in  the  0.4-0. 8  Hz  frequency  band 
for  all  sites  studied.  Similar  correlations  were  noted  for  Novaya  Zemlya  and  Am¬ 
chitka  for  the  0.8-1. 1  Hz,  1.1-1. 5  Hz  and  0. 2-2.0  Hz  bands,  while  Pahute  Mesa 
explosions  exhibit  a  positive  correlation  of  complexity  with  magmtude,depth  in  the 
0.8-1. 1  Hz  band.  Spatial  correlations  of  coda  complexity  were  also  noted  for  the 
Pahute  Mesa  test  site,  with  complexity  increasing  toward  the  center  of  the  test  site. 
No  correlation  with  tectonic  release  F-factor  was  found;  thus,  it  appears  that  after¬ 
shocks  do  not  contribute  significant  high  frequency  energy  to  the  coda,  as  had  been 
suggested  by  Douglas  ( 1984)  on  the  basis  of  analysis  of  the  Pahute  Mesa  explosion 
GREELEY. 

Another  interesting  result  from  the  study  of  Pahute  Mesa  coda  was  the  detection 
of  systematic  azimuthal  patterns  in  the  early  coda  amplitudes  that  mirror  the  patterns 
in  the  direct  arrivals,  albeit  with  a  somewhat  smaller  range  of  variation.  In  a  detailed 
analysis  of  the  direct  signals.  Lynnes  and  Lay  ( 1988)  demonstrated  a  correlation 
between  low  mf  and  fast  travel  times  for  Pahute  Mesa  explosions,  which  they 
attributed  to  strong  focussing/defocussing  effects  in  the  lithosphere  and  in  the  deep 
mantle.  Thus,  it  appears  that  the  near-source  component  of  the  early  coda  is  affected 
by  the  same  deep  mantle  heterogeneity,  rather  than  having  effectively  averaged  out  all 
of  the  propagational  effects.  Similar  correlations  between  coda  and  direct  arrival 
amplitudes  have  been  observed  by  Bullitt  and  Cormier  ( 1984)  and  Balmgardt 
(1985). 

In  this  study,  we  have  digitized  1063  WWSSN  and  Canadian  Seismic  Network 
seismograms  from  32  Yucca  Flat  explosions  to  complement  the  1019  seismograms 
from  25  Pahute  Mesa  tests  analyzed  by  Lay  and  Wf.lc  ( 1987)  and  Lay  ( I987a.b). 
As  in  these  earlier  studies,  we  identify  the  near-source  contribution  to  the  early 
P  coda  by  using  a  large,  well-distributed  station  set.  The  dependence  of  coda 
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complexity  on  various  event  parameters  such  as  depth  and  magnitude  is  determined 
for  the  two  test  sites.  We  also  examine  the  individual  event  average  complexity  values 
in  order  to  look  for  coherent  spatial  variations  within  the  test  sites.  Such  variations 
provide  prime  motivations  for  future  aitempts  to  identify  particular  scattering  struc¬ 
tures.  Finally,  the  azimuthal  patterns  of  the  direct  arrivals  and  coda  for  the  two 
nearby  test  sites  will  be  compared  to  determine  the  relative  importance  of  near-source 
and  deep  mantle  effects. 


Data  Analysis 

The  seismograms  for  Yucca  Flat  events  were  processed  in  a  nearly  identical 
manner  to  those  for  Pahute  Mesa  events,  as  described  in  Lay  and  Welc  ( 1987)  and 
Lay  ( 1987a).  At  least  15  sec  of  noise  and  15  sec  of  signal  were  digitized  and  interpo¬ 
lated  at  a  0.1  sec  sample  interval.  The  seismograms  were  filtered  using  three-pole 
Butterworth  bandpass  filters  with  four  different  passbands:  0. 2-2.0  Hz  (“broad¬ 
band"),  0.4-0.8  Hz  (low-frequency).  0.8—1 . 1  Hz  (mid-frequency),  and  1.1-1. 5  Hz 
(high-frequency).  The  individual  instrument  responses  were  deconvolved  from  the 
low-frequency,  mid-frequency  and  high-frequency  data  sets  before  filtering;  however, 
deconvolutions  of  the  “broadband"  seismograms  were  unstable  at  the  long  period 
end,  so  the  instrument  response  was  not  removed  from  this  data  set.  Since  most  of  the 
stations  are  from  the  WWSSN,  this  results  only  in  a  slight  increase  of  the  variance  in 
the  broadband  results. 

In  order  to  characterize  the  complexity  of  the  coda,  the  normalized  energy  tempo¬ 
ral  centroid  (McLaughlin  and  Anderson,  1985)  of  the  first  15  sec  of  the  squared 
analytic  envelope  (Figure  I)  was  computed.  The  energy  temporal  centroid  (Chp)  is 
defined  for  each  bandpass  filtered  (bp)  seismogram  by: 

"15 

Chn  =  I  s'(t)t  dt.  s:(t )  dt. 

Jo  Jo 

where  s'~(t)  is  the  square  of  the  filtered  signal  envelope.  Larger  centroid  values  indi¬ 
cate  a  higher  proportion  of  energy  later  in  the  signal.  A  correction  for  noise  level  was 
applied  by  subtracting  the  root-mean-square  (RMS)  of  the  noise  window  from  the 
envelope  before  computing  Chp;  this  generally  reduced  all  the  event  average  centroids 
by  comparable  amounts.  Station/path  corrections  were  determined  by  averaging  the 
centroids  for  all  events  from  a  given  test  site  recorded  at  a  given  station.  When 
computing  an  event  average  value  of  C^,,.  the  station  corrections  were  applied,  and 
the  corrected  centroids  were  azimuthally  averaged  in  7.5  windows  to  reduce  bias 
from  uneven  azimuthal  distribution.  The  event  centroids  for  each  test  site  were  then 
demeaned,  so  that  the  values  reported  below  arc  deviations  from  the  average. 

The  RMS  of  the  direct  arrival  (0  5  sec)  and  the  coda  (5-15  sec)  were  also  com¬ 
puted  for  each  seismogram.  The  logarithm  of  the  ratio  RMS'  ”  RMS0^ '  provides 
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Figure  1 

Example  of  the  data  analysis  procedure.  The  top  trace  is  the  Pwavetram  with  its  analytic  envelope,  filtered 
from  0.2  -2.0  Hz  for  the  Yucca  Flat  event  TAN  at  station  ATL.  RMS  amplitudes  are  calculated  for  the  first 
5  sec  (direct  arrival)  and  the  5  15  sec  window  (early  coda).  The  bottom  trace  shows  the  squared  envelope, 
before  the  noise  correction  (solid  line)  and  after  the  noise  correction  (dashed  line),  and  the  effect  of  the 
noise  correction  on  the  energy  temporal  centroid 


another  measure  of  the  complexity  of  the  coda,  with  higher  values  indicating  rela¬ 
tively  enhanced  coda  levels.  Event  averages  for  log(RMS5  ^  'RMS0  5)  were  calcu¬ 
lated  in  the  same  manner  as  for  the  centroids.  The  log(RMS5l5/RMSa  5)  measure 
was  used  to  compare  with  the  behavior  of  the  direct  arrival  (log(RMS°~5))  and  early 
coda  (log(RMS5  ,5));  however  the  centroids  will  be  preferred  in  this  study  when 
considering  event  average  complexity,  since  it  proved  difficult  to  apply  a  stable  noise 
correction  to  the  log(  RMS5  ''/RMS^5)  measurements.  The  correlation  coefficients. 
r.  of  the  event  average  centroids  (without  noise  correction)  vs.  log(RMS5  l5/RMS°~5) 
are  r  =0.94  in  the  low-frequency  band,  r  =  0.84  in  the  mid-frequency  band,  and 
r  =  0.94  in  the  high-frequency  band  for  Yucca  Flat  events.  The  broadband  correla¬ 
tion  coefficient  is  only  r  =  0.72,  which  may  be  due  to  the  fact  that  the  seismograms 
were  not  deconvolved  in  this  case,  fn  general,  however,  these  two  quite  different  types 
of  coda  complexity  measures  gave  very  similar  results,  as  was  the  case  for  the  Pahute 
Mesa  date  (Lay  and  Welc.  1987;  Lav.  1987a). 

Event  Average  Complexity 

Correlation  with  mb  and  depth.  The  event  average  centroids  in  each  frequency 
band  for  Yucca  Flat  were  correlated  wiih  magnitude,  burial  depth  and  depth  below 
the  water  table.  However,  the  only  parameter  that  gave  a  correlation  coefficient 
above  0.5  in  any  of  the  frequency  bands  was  magnitude.  The  correlations  of  coda 
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complexity  with  magnitude  and  with  burial  depth  in  each  frequency  hand  for  both 
Yucca  Flat  and  Pahute  Mesa  events  are  listed  in  Table  I.  Variance-weighted  lincai 
regressions  of  coda  complexity  vs.  magnitude  and  depth  for  both  test  sites  are  shown 
in  Figure  2  for  the  low-frequency  -and  mid-frequency  bands. 

both  test  sites  have  little  magnitude  dependence  in  the  broadband  and  high- 
frequency  band  cases.  In  the  mid-frequency  band,  however,  the  Pahute  Mesa  events 
show  a  positive  correlation  of  complexity  with  magnitude  (r  =  0.62)  and  depth 
(r  =  0.71).  There  is  no  correlation  with  magnitude  for  Yucca  Flat,  but  a  similar  trend 
with  burial  depth,  although  several  outliers  lead  to  a  poor  correlation  coefficient 
(/■  =  0.27).  These  correlations  may  result  from  variation  of  the  spectrum  of  the  direct 
arrival  in  this  band  due  to  pP  interference.  For  larger,  more  deeply  buried  events,  the 
increase  in  pP  delay  time  would  shift  the  spectral  scalloping  from  pP  interference 
toward  the  middle  of  the  passband.  This  effect  can  be  predicted  b>  computing  syn¬ 
thetic  spectra  for  a  large,  deep  event  and  a  small,  shallow  event  Two  Pahute  Mesa 
explosions  were  used  for  these  computations  (Figure  3).  because  the  teleseismic  P 
waveforms  for  tests  at  this  site  have  been  extensively  studied  (Las.  1985).  The 
spectrum  of  the  coda  should  not  have  such  pronounced  scalloping  due  to  pP  (Glpta 
and  Blandford.  1987)  and  thus  the  coda  amplitudes  in  this  passband  would  appear 
relatively  enriched. 

In  order  to  demonstrate  this  with  actual  data,  we  compared  the  spectral  ampli¬ 
tudes  of  a  15s  noise  window  to  two  1 5  s  windows  in  the  P  wavetrain.  one  including 
the  direct  arrival  (0-12.8  s)  and  one  without  the  direct  arrival  (5  1 7.8  s).  This  method 
of  presentation  ( Figure  4)  mitigates  the  possible  effects  of  noise  bias  between  a  large 
and  small  explosion.  COMMODORE  is  a  large,  deeply  buried  Yucca  Flat  test,  while 
PIRANHA  is  a  relatively  small,  shallow  event,  with  an  average  centroid  that  is  0.7  s 
less  than  that  of  COMMODORE  in  the  mid-licquency  band  The  spectral  null 
related  to  pP  interference  is  clearly  visible  in  the  0  1 2.8  s  window,  which  is  dominated 
by  the  direct  phase,  for  both  events  As  in  the  s\nthenc  case  (Figure  3)  the  null  is 
shifted  into  the  mid-frequency  band  for  the  deeper  test.  This  depletes  the  direct 
arrival  relative  to  the  coda,  in  which  the  spectral  null  ts  much  subdued  as  seen  in  the 
pure  coda  spectra  (5.0-17.8  s).  Nevertheless,  the  null  does  exert  a  small  effect  on  the 
spectral  amplitudes,  suggesting  that  the  coda  has  some  component  of  P  P  scattering. 


Tabic  I 


Correlation  coefficients  of  C\r 

with  magnnutit’ 

and  depth  of  hunal 

Test  Site 

Parameter 

0  4  UK  11/ 

08  111!/ 

II  ISM- 

0  2  2  0  It/ 

Yucca  Flat 

-061 

006 

-0  19 

0  3S 

Pahute  Mesa 

-0  94 

0.62 

-0  27 

-  0  34 

Yucca  Flat 

depth 

-0  33 

0.27 

0  09 

-0  06 

Pahute  Mesa 

depth 

-0  80 

071 

-0  38 

-000 
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Yucca  Flat  Pahute  Mesa  Yucca  Flat  Pahute  Mesa 


Figure  2 

Correlations  of  relative  event  average  centroids  (AC*,,)  with  magnitude  and  burial  depth  in  the  low- 
frequeno  and  mid-frequency  bands  for  Yucca  Flat  events  and  Pahute  Mesa  events 


log(frequency) 

Figure  3 

Synthetic  amplitude  spectra,  calculated  using  r*  =  0  75  sec  and  a  short-period  WWSSN  instrument  for  a 
large,  deeply  buned  event  (BOXCAR  yield  =  1300  kt.  depth  =  I  158  km)  and  a  smaller,  less  deeply 
buried  event  ( HALFBEAK;  yield  =  365  kt.  depth  =  0.819  km)  at  Pahute  Mesa  The  increase  in  pP  delay 
time  for  the  large  event  shifts  the  scalloping  effect  of  pP  interference  into  the  mid-frequency  band 
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Figure  4 

Ratio  of  seismogram  spectral  amplitude  to  noise  for  a  large,  deep  event  (COMMODORE)  and  a  small, 
shallow  event  (PIRANHA)  at  Yucca  Flat  The  dashed  line  indicates  the  0-12  8  s  window,  which  is 
dominated  by  the  direct  phase,  and  the  5-17.8  :>  window,  which  does  not  include  the  direct  phase 


Another  possibility  is  frequency  dependent  defocussing  of  the  direct  arrival  for 
events  in  the  center  of  the  mesa,  where  most  of  the  larger  events  were  located.  This 
second  possibility  will  be  tested  later  by  comparing  spatial  variations  of  complexity 
and  azimuthal  patterns  in  the  different  frequency  bands. 

In  the  low-frequency  band,  both  test  sites  show  a  negative  correlation  with  mag¬ 
nitude  (r  =  -0.94  for  Pahute  Mesa  and  r  =  —  0.61  for  Yucca  Flat),  with  smaller 
magnitude  events  having  relatively  high  coda  levels.  For  Pahute  Mesa,  the  low- 
frequency  band  has  a  similar  dependence  on  depth  of  burial  (r  =  -0.80).  For  Yucca 
Flat  events,  the  dependence  of  complexity  on  depth  is  fairly  weak  (r  =  -0.33), 
suggesting  that  the  correlation  of  complexity  with  magnitude  cannot  be  ascribed  to 
depth  of  burial  effects  alone,  but  it  is  important  to  note  that  the  Pahute  Mesa  events 
span  a  much  greater  range  in  burial  depth. 

One  possibility  for  these  correlations  is  the  effect  of  pP  interference.  However,  the 
parallel  spectral  curves  in  the  0.4-0. 8  Hz  band  for  both  explosions  in  Figure  4  indicate 
that  this  band  is  relatively  free  of  such  interference.  Furthermore,  the  effect  of  inter¬ 
ference  should  be  to  decrease  the  direct  arrivals  of  the  larger  explosions,  as  in  the 
mid-frequency  band,  which  would  result  in  a  trend  opposite  of  that  which  is  observed. 
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Another  possibility  is  that  the  negative  correlation  of  complexity  with  size  is  due 
to  greater  noise  contamination  in  this  band  for  smaller  explosions.  The  signal/noise 
ratio  is  indeed  smaller  for  PIRANHA  than  for  COMMODORE  (Figure  4),  but  it  is 
still  on  average  fairly  high  within  this  band  for  both  tests,  dipping  below  2  only  for 
the  coda  of  PIRANHA.  The  effect  of  noise  is  further  reduced  by  application  of  the 
noise  correction  to  the  centroids.  Also,  although  noise  will  in  general  increase  the 
centroids,  noise  that  is  in  phase  with  the  direct  arrival  decreases  the  centroids,  further 
mitigating  the  overall  effects  of  noise.  These  effects  are  also  reduced  by  the  inclusion 
of  low-noise  stations,  for  which  the  effects  of  noise  are  negligible.  Thus,  it  seems 
unlikely  that  noise  contamination  alone  could  account  for  the  substantial  variation 
in  centroid  with  magnitude  in  the  low-frequency  band 

Spatial  variations  of  complexity.  The  spatial  variation  of  event  average  complexity 
in  all  four  passbands  for  Yucca  Flat  is  shown  in  Figure  5.  The  coherent  variations  in 
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Figure  5 

Spatial  variation  ot  event  average  complexity  (ACV)  for  the  Yucca  Flat  test  site  in  the  four  frequency 
bands  Diamonds  indicate  smaller  than  average  coda  complexity,  and  crosses  larger  than  average  coda 

complexity 
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complexity  with  location  are  quite  striking.  The  least  coherent  variations  occu;  m  the 
low-frequency  band;  thus,  the  complexity  in  this  band  appears  to  be  influenced  more 
by  magnitude  than  by  the  source  environment.  The  most  coherent  variations  occur 
in  the  mid-frequency  band.  This  result  is  consistent  with  the  frequency  dependence  of 
spatial  variations  for  Pahute  Mesa,  where  the  mid-frequency  band  had  the  highest 
correlation  with  distance  from  the  center  of  the  mesa  (r  =  -0.72  vs.  r  =0.24  for  the 
low-frequency  band  and  r  =  —  0.16  for  the  high-frequency  band).  The  actual  spatial 
patterns  are  also  somewhat  different  among  the  different  frequency  bands:  events  in 
the  northwest  part  of  the  test  array  have  very  low  complexity  in  the  mid-frequency 
band,  but  relatively  high  complexity  in  the  high-frequency  band;  events  in  the  north- 
central  part,  on  the  other  hand,  have  high  complexity  in  the  mid-frequency  band,  but 
average  complexity  in  the  high-frequency  band.  The  broadband  map  resembles  a 
subdued  version  of  the  mid-frequency  band  map.  due  to  the  predominance  of  mid¬ 
frequency  band  energy  in  the  undeconvolved  broadband  signals  and  the  compen¬ 
sating  trends  in  the  different  frequency  bands.  The  difference  in  patterns  among 
frequency  bands  suggests  a  frequency  dependence  in  focussing  of  the  direct  arrival  or 
in  scattering  contributions  to  the  coda. 

The  spatial  variations  of  coda  complexity  imply  a  correlation  with  some  properly 
of  the  source  region.  Ferguson  s  (1982)  model  for  Yucca  Flat  basin,  based  on 
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borehole  data  and  gravity  modeling,  is  shown  in  Figure  6  with  the  event  average 
complexity  values  for  the  mid-frequency  band  superimposed.  The  contours  indicate 
the  depth  to  the  Tertiary-Paleozoic  contact.  The  events  with  the  highest  complexity 
lie  above  the  basin  axis  at  the  north  end,  just  to  the  west  of  the  Yucca  Fault.  The 
events  with  the  simplest  coda  levels  are  situated  above  the  east  flank  of  a  horst  in  the 
basement.  This  suggests  that  either  the  basin  shape  or  the  thickness  of  Tertiary  tuffs 
controls  coda  complexity  in  this  passband. 


Azimuthal  Patterns  of  Amplitudes  and  Complexity 

Slowly  varying  azimuthal  patterns  of  mf'  and  travel  time  residuals  have  been 
observed  for  events  at  both  the  Pahute  Mesa  and  Yucca  Flat  test  sites  (Lynnes  and 
Lay,  1988).  Correlations  between  the  travel  time  and  amplitude  variations  indicate 
that  a  significant  component  is  due  to  focussing  and  defocussing  by  velocity  hetero¬ 
geneities  in  the  mantle.  Minster  et  al.  ( 1981)  and  Taylor  ( 1983)  found  localized 
high  velocity  regions  in  the  upper  mantle  to  the  northeast  of  Pahute  Mesa,  which  can 
produce  defocussing  that  accounts  for  much  of  the  difference  between  the  Pahute 
Mesa  and  Yucca  Flat  amplitude  patterns  (  Lynnes  and  Lay,  1988;  Cormier.  1987). 
However,  the  patterns  for  the  two  test  sites  share  a  large  common  component,  which 
Lynnes  and  Lay  ( 1988)  inferred  to  be  due  to  defocussing  and  focussing  by  deep 
<  >  200  km)  heterogeneity.  !t  is  thus  of  interest  to  determine  the  relative  influence  of 
such  heterogeneities  on  the  direct  arrival  and  the  early  coda. 

In  order  to  determine  the  average  azimuthal  patterns  for  the  direct  arrivals 
(log(RMS“  ')).  early  coda  (log(RMS'  ls))  and  complexity  (log(RMS5  ls  RMSl>  5)), 
the  station  averages  for  each  site  were  calculated  using  a  least-squares  technique  that 
inverts  for  the  best  relative  source  sizes  to  simultaneously  minimize  the  variance 
about  the  station  means  (  Larry  Ruff,  personal  communication).  In  order  to  bring  out 
the  slowly  varying  component  of  the  azimuthal  patterns,  curves  of  the  form 
A  sin[2((f  -  ft,,)]  were  regressed  through  the  resultant  station  averages  as  a  function 
of  azimuth  (0).  The  particular  form  of  the  regression  curve  is  not  physically  moti¬ 
vated,  but  is  intended  to  parameterize  the  long-wavelength  variations  to  facilitate 
comparisons  between  events  and  test  sites.  The  F-ratio  is  also  calculated  as  described 
in  Lay  et  al.  ( 1984)  to  provide  an  indicator  of  the  coherence  and  significance  of  the 
long-wavelength  component.  The  results  for  the  broadband  set  are  shown  for  Pahute 
Mesa  and  Yucca  Flat  in  Figure  7. 

The  direct  arrivals  (log(  RMS'* '))  show  azimuthal  patterns  for  both  Pahute  Mesa 
and  Yucca  Flat,  with  a  somewhat  stronger  and  more  coherent  pattern  for  Pahute 
Mesa.  The  correlation  in  the  log(RMSu  5)  station  averages  between  the  test  sites  is 
r  =  0.83,  while  it  is  r  =  0.73  for  the  log(  RMS'  1 ')  values.  Thus,  a  significant  compo¬ 
nent  of  the  variation  is  common  to  the  two  test  sites.  Lynnes  and  Lay  (1988) 
obtained  a  similar  result  for  mJhh  and  interpreted  the  higher  amplitude  of  variation  for 
Pahute  Mesa  as  a  result  of  focussing  and  defocussing  by  lateral  velocity  gradients  in 
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I  igure  7 

Sin  20  regressions  for  Yucca  Flat  and  Pahule  Mesa  station  jseraged  logi  RMS"  'i  logi  RMS'  'i  and 
log(  RMS5  ”  RMS"”')  in  the  0  2  2  0  Hz  passband  The  sine  curses  are  dashed  vs  here  not  statistical 

significant  at  the  95°t»  level 

the  crust  and  shallow  mantle  under  the  mesa.  The  log(  RMS"  )  pattern  for  Pahute 
Mesa  is  stronger  than  the  iog(  RMS'  ")  pattern  lor  the  same  site,  indicating  that  the 
focussing  defocussing  in  the  Pahute  Mesa  direct  signals  is  less  pronounced  for  the  earlv 
coda  due  to  the  averaging  characteristics  of  the  scattered  wavefield  The  log(  RMS'  1 ') 
pattern  for  Yucca  Flat,  however,  has  a  stronger  long-wavelength  component 
(F-ratio  =  7.6)  with  a  slightly  larger  range  than  the  pattern  for  the  direct  arrivals 
(F-ratio  =  3.9)  from  that  site,  although  the  orientation  of  the  patterns  are  quite  simi¬ 
lar.  Th  s  behavior  is  present  in  many  sine  curve  regressions  of  logi  RMS"  ')  and 
log(  RMS'  ’')  of  the  individual  events  as  well.  Thus,  it  seems  that  both  the  direct  arri¬ 
vals  am:  the  coda  arc  affected,  more  or  less  equally,  by  deep  mantle  heterogeneity,  while 
the  decreased  significance  of  the  direct  arrival  sine  curves  appears  to  be  due  to  increased 
scatter  in  the  associated  station  averages.  This  might  result  from  very-near-source  or 
near-receiver  heterogenc.tv.  The  Yucca  Flat  events  were  detonated  close  to  large 
basement  offsets  on  the  Yucca  Fault  ( Figure  6).  which  may  cause  triplications  and 
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caustics  affecting  the  direct  arrivals  and  early  coda  (Stead  and  Helmberger,  1986), 
In  any  case,  the  homogenization  of  the  early  coda  seems  to  have  little  effect  on  sub¬ 
duing  the  predominant  amplitude  variations  induced  by  earth  structure  at  both  sites. 

Sine  curves  were  also  regressed  through  the  inverse  complexity  measure 
log(RMS°~5/RMS5~15)  for  Pahute  Mesa  and  Yucca  Flat  (Figure  7).  As  might  be  ex¬ 
pected,  the  Pahute  Mesa  sine  curve  is  significant  (F-ratio  =  5.2),  due  to  the  larger  range 
of  variations  for  log(RMS°~5)  than  for  (RMS515).  The  Yucca  Flat  pattern,  on  the 
other  hand,  does  not  have  any  long-wavelength  variation,  indicating  that  the  direct 
arrival  and  the  phases  making  up  the  coda  are  affected  similarly  by  deep  mantle 
structure. 
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Figure  8 

Map  of  Yucca  Flat  showing  the  /-  -ratios  and  azimuths  of  the  negative  lobes  of  sin  20  regressions  for 
individual  events  in  the  low-frequency,  mid-frequency  and  high-frequency  bands.  The  first  column  shows 
results  for  loglRMS®  ').  the  second  shows  logtRMS'  ")  The  length  of  the  arrow  is  proportional  to  the 
/'-ratio  Arrows  for  events  with  an  /'-ratio  lower  than  3  arc  not  shown 
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In  order  to  further  explore  the  frequency  dependence  of  the  focussing  and  or 
scattering  effects,  sine  curves  were  regressed  through  log(  RMSU  5).  log(  RMS5  l5)  and 
logfRMS^’/RMS5  l5)  measures  of  individual  events  from  both  test  sites  in  all  four 
frequency  bands.  The  F-ratios  and  the  directions  of  the  negative  lobes  of  the  direct 
arrival  and  early  coda  patterns  in  the  low-frequency,  mid-frequency  and  high-fre¬ 
quency  bands  for  Yucca  Flat  are  shown  in  Figure  8.  As  with  the  broadband  results, 
the  regressions  of  log(  RMS1*"')  values  have  systematically  smaller  F-ratios  than  the 
!og(  RMS5  l5)  measurements:  however,  the  amplitudes  and  azimuths  are  quite  similar 
for  a  given  frequency  band.  In  addition,  there  is  a  systematic  clockwise  rotation  in  the 
long-wavelength  patterns  with  frequency  from  2  4  (azimuth  of  the  negative  lobe) 
in  the  low-frequency  band  to  11  -  13  in  the  mid-frequency  band  to  17  -20  in  the 
high-frequency  band. 

The  /"-ratios  and  azimuths  of  individual  event  patterns  in  the  inverse  complexity 
for  all  four  frequency  bands  are  shown  in  Figure  9  for  both  test  sites.  There  were  too 
few  significant  sine  curves  to  deduce  any  frequency  dependence  in  azimuthal  patterns 
of  complexity  for  Yucca  Flat,  but  a  spatial  dependence  does  appear  in  that  those  few 
sine  curves  that  are  significant  tend  to  be  for  events  at  the  southern  end  of  the  test 
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Figure  9 

Map  of  the  Nevada  Tesl  Site  showing  the  f-ratios  and  azimuths  of  the  negative  lobes  of  sin  20  regressions 
through  loglRMS^'  RMS'  "l  (inverse  complexity)  for  individual  events  in  the  broadband,  low- 
frei  :ency,  mid-frequency  and  high-frequency  ranges 
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array  On  the  other  hand.  Pahute  Mesa,  events  clearly  exhibit  a  frequency  depen¬ 
dence.  as  there  are  more  events  with  coherent  patterns  tn  the  mid-frequency  band 
than  in  any  other  band.  Thus,  it  appears  that  it  is  the  signal  in  this  passband  that 
produces  the  significant  azimuthal  patterns  for  Pahute  Mesa  complexity  observed  by 
Lay  and  Wei.c  ( 1987)  for  events  in  the  middle  of  the  mesa.  Furthermore,  almost  all 
of  the  events  with  significant  patterns  in  this  passband  are  in  the  central  part  of  the 
mesa,  supporting  the  hypothesis  that  frequency  dependence  of  focussing/defocussing 
or  scattering  of  the  direct  arrival  plays  a  role  in  producing  these  patterns.  Because  the 
large  events  tend  to  be  located  in  the  center  of  the  mesa,  frequency-dependent 
focussing  may  also  help  to  explain  the  positive  correlation  of  complexity  with  magni¬ 
tude  in  the  mid-frequency  band,  rather  than  simply  pP  variations. 


Discussion 

One  of  the  most  difficult  observations  to  explain  from  this  paper  and  Lav  (  1987a) 
is  the  dependence  of  complexity  on  magnitude  All  sites  studied  so  far  have  exhibited 
enriched  long  period  (0.4  4)  8  Hz)  energy  in  the  coda,  relative  to  the  direct  arrival  for 
smaller  events.  At  Pahute  Mesa,  this  effect  could  in  fact  be  a  dependence  on  depth  of 
burial.  Lay  ( 1987a)  suggests  that  a  negative  correlation  of  long  period  complexity 
with  depth  could  ensue  from  increased  Rayleigh-to-P  scattering  for  the  shallower 
events.  However,  the  depth  dependence  is  slightly  weaker  than  the  magnitude  depen¬ 
dence  at  Pahute  Mesa,  and  is  quite  a  bit  weaker  than  the  magnitude  dependence  at 
Yucca  Flat,  suggesting  that  some  variable  more  directly  related  to  source  size  con¬ 
trols  long-period  complexity.  Such  an  effect  would  seem  to  be  intrinsically  nonlinear, 
and  as  vet  there  does  not  seem  to  be  an  adequate  explanation  for  this  phenomenon. 
One  possibility  is  that  the  cavity  size  influences  the  excitation  of  near-field  Rayleigh 
waves,  which  scatter  into  the  teleseismic  wavefield,  or  that  asymmetry  exists  in  the 
spectrum  of  energy  radiated  horizontally  relative  to  that  radiated  vertically.  Numer¬ 
ical  simulations  are  needed  to  explore  these  possibilities. 

The  spatial  variations  in  complexity  found  for  both  Pahute  Mesa  and  Yucca  Flat 
strongly  suggest  that  information  about  the  source  region  structure  is  contained  in 
the  teleseismic  P  wave  coda.  These  variations  are  slightly  dependent  on  the  frequency 
band,  being  particularly  coherent  in  the  mid-frequency  band  at  both  test  sites.  The 
spatial  dependence  of  event  complexity  at  Pahute  Mesa  is  difficult  to  separate  from 
the  size  dependence  since  the  largest  events  are  located  near  the  center  of  the  mesa. 
However,  the  strong  azimuthal  patterns  of  complexity  and  log(RMS°  s)  for  events 
near  the  center  of  the  mesa  suggest  that  defocussing  of  the  direct  arrivals  is  part 
of  the  cause  of  the  enhanced  complexity  of  the  coda  in  the  mid-frequency  band. 
The  spatial  dependence  in  the  mid-frequency  band  for  Yucca  Flat  is  slightly  easier 
to  interpret  since  event  complexity  in  this  band  is  not  complicated  by  depth  or 
magnitude  dependence.  Clearly,  some  near-source  factor  affects  complexity.  Possible 
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candidates  are  the  thickness  of  the  Tertiary  tuffs,  which  may  give  rise  to  larger 
Rayleigh  waves  in  the  center  of  the  basin  (  Figure  6).  or  the  basin  shape  itself,  which 
may  defocus  the  direct  arrivals  to  the  north  for  events  in  the  center  of  the  basin. 
However,  it  is  difficult  to  determine  what  effect  these  two  factors  would  have  without 
sophisticated  finite-difference  modeling  (e  g..  Stead  and  Hllmberger.  1986). 

While  the  event  average  complexity  results  point  to  near-source  information  in 
the  ratio  of  the  coda  to  the  direct  arrival,  the  azimuthal  patterns  of  the  early  coda  for 
Yucca  Flat  emphasize  the  sizable  contribution  of  a  deep  mantle  effect.  The  most 
surprising  result  is  the  high  significance  of  the  long  wavelength  component  for  the 
coda  relative  to  the  direct  arrival.  This  behavior  contrast';  with  that  for  Pahute  Mesa, 
where  the  direct  arrivals  seem  to  be  more  strongly  defocussed  to  the  north  by  rela¬ 
tively  shallow  structure  than  the  coda.  The  increased  scatter  in  the  Yucca  Flat  direct 
signals  suggests  that  the  direct  arrivals  from  the  individual  events  pass  through  very 
shallow  structure  that  has  little  influence  on  the  scattered  phases  making  up  the  coda 
(Figure  10).  Thus,  the  coda  comprises  arrivals  that  leave  the  source  region  with  a 


Pahute  Yucca 

Mesa  Flat 


Schematic  diagram  shnw.nt.’  the  effects  of  near-source  and  deep  mantle  heterogeneity  at  the  Nevada  Test 
Site,  Direct  arrivals  from  Pahute  Mesa  are  defocussed  by  a  high-velocity  structure  that  does  not  affect  the 
scattered  phase.  Direct  arriv.i  ,  from  Yucca  Flat  arc  complicated  by  triplications  and  caustics  resulting 
from  the  basin  structure  All  of  the  arrivals  arc  then  affected  similarly  by  a  deep,  regional  mantle  hetero¬ 
geneity.  which  imparts  a  common  azimuthal  pattern  to  direct  arrivals  and  coda  from  both  test  sites 
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smaller  amplitude  variability  than  the  direct  downgoing  energy.  Then  the  entire 
signal  encounters  the  deep  mantle  heterogeneity  inferred  by  Lynnes  and  Lay  ( 1988), 
which  imparts  a  common  azimuthal  pattern.  The  direct  arrivals  have  the  same  super¬ 
imposed  variation,  but  the  significance  of  the  long  wavelength  component  is  reduced 
by  the  additional  scatter  from  very  near-source  or  near-receiver  heterogeneity. 

The  presence  of  strong  azimuthal  patterns  in  the  early  P  coda  for  Yucca  Flat 
events  raises  the  question  of  whether  the  early  P  coda  is  a  more  precise  estimator  of 
source  size  than  the  direct  arrival  for  Yucca  Flat  events.  To  answer  this  question,  the 
event  standard  deviations  were  averaged  for  the  three  measures  mf.  log(RMS°“5), 
and  log(RMS5'15):  the  average  standard  deviations  are  0.14  +  .02  for  0.15  ±  .02 
for  log(RMS"  '),  and  0.18  ±  .02  for  log(  RMS'  l5).  Thus,  the  early  P  coda  actually 
has  somewhat  higher  variance  than  the  direct  arrival  for  Yucca  Flat  events.  While 
this  difference  may  be  partially  due  to  the  increased  effect  of  noise,  most  of  the 
variance  in  the  early  coda  would  seem  to  result  from  the  long  wavelength  component 
produced  by  deep  mantle  heterogeneity. 

Although  the  network-averaged  coda  measurements  from  Yucca  Flat  events  are 
less  precise  than  the  mh  measurements.  Baumgardt  ( 1985)  found  that  the  standard 
deviation  of  P  coda  across  NORSAR  for  Semipalitinsk  events  was  even  less  than  that 
of  the  noise  across  the  array.  This  suggests  that  array-averaged  coda  measurements 
are  more  precise  than  network-averaged  P  or  P  coda.  Indeed,  azimuthal  patterns  in 
the  coda  such  as  those  from  NTS  would  produce  a  large  degree  of  scatter  in  network- 
averaged  measurements,  but  very  little  additional  scatter  at  a  single  array,  leading  to 
a  comparable  or  higher  precision  at  the  latter. 

Even  more  important  than  the  precision  of  the  different  measures,  however,  is 
their  accuracy.  Clearly,  measurements  of  the  coda  from  NTS  at  a  single  station  or 
array  have  a  high  probability  of  falling  on  one  of  the  lobes  of  the  azimuthal  pattern 
from  NTS.  (NORSAR.  for  instance,  occurs  near  the  minimum  of  the  Pahute  Mesa 
and  Yucca  Flat  azimuthal  amplitude  patterns.)  Consequently,  single  arrays  may 
record  coda  magnitudes  with  very  high  precision  by  averaging  out  the  near-receiver 
effects,  but  cannot  average  out  effects  imparted  by  mantle  heterogeneity  near  the 
source  region,  resulting  in  an  absolute  bias. 


Conclusions 

Data-intensive  analysis  of  teleseismic  P  wave  signals  has  been  used  to  isolate 
near-source  contributions  to  the  early  coda  from  NTS  events.  The  0.4--0.8  Hz  fre¬ 
quency  band  has  a  strong  negative  correlation  between  coda  complexity  and  magni¬ 
tude  at  both  the  Pahute  Mesa  and  Yucca  Flat  test  sites,  with  smaller  events  having 
higher  coda  levels.  Effects  related  to  depth  of  burial  alone  do  not  seem  sufficient  to 
explain  this  correlation,  since  the  correlation  of  complexity  with  depth  for  Y ucca  Flat 
is  fairly  low. 
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smaller  amplitude  variability  than  the  direct  downgoing  energy.  Then  the  entire 
signal  encounters  the  deep  mantle  heterogeneity  inferred  by  Lynnes  and  Lay  ( 1988), 
which  imparts  a  common  azimuthal  pattern.  The  direct  arrivals  have  the  same  super¬ 
imposed  variation,  but  the  significance  of  the  long  wavelength  component  is  reduced 
by  the  additional  scatter  from  very  near-source  or  near-receiver  heterogeneity. 

The  presence  of  strong  azimuthal  patterns  in  the  early  P  coda  for  Yucca  Flat 
events  raises  the  question  of  whether  the  early  P  coda  is  a  more  precise  estimator  of 
source  size  than  the  direct  arrival  for  Yucca  Flat  events.  To  answer  this  question,  the 
event  standard  deviations  were  averaged  for  the  three  measures  mf,  log(RMSa5), 
and  log(  RMS5  IS):  the  average  standard  deviations  are  0.14  +  .02  for  mahb ,  0.15  ±  .02 
for  log(RMS"  '),  and  0.18  ±  .02  for  logfRMS'  l5).  Thus,  the  early  P  coda  actually 
has  somewhat  higher  variance  than  the  direct  arrival  for  Yucca  Flat  events.  While 
this  difference  may  be  partially  due  to  the  increased  effect  of  noise,  most  of  the 
variance  in  the  early  coda  would  seem  to  result  from  the  long  wavelength  component 
produced  by  deep  mantle  heterogeneity. 

Although  the  network-averaged  coda  measurements  from  Yucca  Flat  events  are 
less  precise  than  the  mh  measurements.  Baumgardt  ( 1985)  found  that  the  standard 
deviation  of  Pcoda  across  NORSAR  for  Semipalitinsk  events  was  even  less  than  that 
of  the  noise  across  the  array.  This  suggests  that  array-averaged  coda  measurements 
are  more  precise  than  network-averaged  P  or  P  coda.  Indeed,  azimuthal  patterns  in 
the  coda  such  as  those  from  NTS  would  produce  a  large  degree  of  scatter  in  network- 
averaged  measurements,  but  very  little  additional  scatter  at  a  single  array,  leading  to 
a  comparable  or  higher  precision  at  the  latter. 

Even  more  important  than  the  precision  of  the  different  measures,  however,  is 
their  accuracy.  Clearly,  measurements  of  the  coda  from  NTS  at  a  single  station  or 
array  have  a  high  probability  of  falling  on  one  of  the  lobes  of  the  azimuthal  pattern 
from  NTS.  (NORSAR.  for  instance,  occurs  near  the  minimum  of  the  Pahute  Mesa 
and  Yucca  Flat  azimuthal  amplitude  patterns.)  Consequently,  single  arrays  may 
record  coda  magnitudes  with  very  high  precision  by  averaging  out  the  near-receiver 
effects,  but  cannot  average  out  effects  imparted  by  mantle  heterogeneity  near  the 
source  region,  resulting  in  an  absolute  bias. 


Conclusions 

Data-intensive  analysis  of  teleseismic  P  wave  signals  has  been  used  to  isolate 
near-source  contributions  to  the  early  coda  from  NTS  events.  The  0.4-0. 8  Hz  fre¬ 
quency  band  has  a  strong  negative  correlation  between  coda  complexity  and  magni¬ 
tude  at  both  the  Pahute  Mesa  and  Yucca  Flat  test  sites,  with  smaller  events  having 
higher  coda  levels.  Effects  related  to  depth  of  burial  alone  do  not  seem  sufficient  to 
explain  this  correlation,  since  the  correlation  of  complexity  with  depth  for  Y ucca  Flat 
is  fairly  low. 
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Inversion  of  P  Coda  for  Isotropic  Scatterers  at  the  Yucca  Flat  Test  Site 
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Abstract 

Although  many  studies  have  addressed  the  statistical  properties  of  P  coda,  there 
have  been  few  attempts  to  analyze  coda  in  a  deterministic  framework.  In  this  analysis, 
1050  teleseismic  P  wave  seismograms  from  32  Yucca  Flat  explosions  are  inverted  for 
the  locations  of  any  isotropic  point  scatterers  in  the  source  region  that  contribute 
coherent  arrivals  to  the  coda.  A  rectangular  grid  of  potential  scattering  locations  is 
constructed  around  the  source  array  and  the  semblance  (a  measure  of  waveform 
ensemble  coherence)  for  each  point  is  computed  by  slant-stacking  the  full  set  of 
seismograms  at  various  moveout  velocities.  The  potential  "resolution",  or  sensitivity, 
of  the  method  for  this  data  set  is  tested  by  constructing  a  corresponding  suite  of  syn¬ 
thetic  seismograms  with  direct  arrivals  plus  arrivals  from  various  point  scatterers  and 
calculating  the  semblance  at  each  point  on  the  grid.  The  particular  source-station 
geometry  produces  better  sensitivity  near  the  array.  Sensitivity  is  enhanced  parallel  to 
the  trend  of  the  elongated  source  array  and  diminished  in  the  orthogonal  direction. 
Stacking  the  data  for  a  velocity  of  2.5  km/s  yields  relatively  high  semblance  values  in 
die  central  portion  as  well  as  about  10  km  to  the  ENE,  and  moderately  high  semblance 
values  about  10  km  to  the  WNW.  The  latter  two  regions  are  associated  with  topo¬ 
graphic  highs.  Simulations  using  synthetics  with  random  arrivals  in  the  coda  indicate 
that  incomplete  suppression  of  random  noise  may  explain  the  high  semblance  values  in 
the  central  region,  but  probably  not  in  the  outlying  regions.  However,  simulations  of  a 
circular  distribution  of  random  scatterers  within  10  km  of  the  test  site  produce  high 
semblance  about  10  km  to  the  ENE,  similar  to  that  seen  for  the  data,  due  to  variable 
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sensitivity.  Thus,  two  possible  explanations  for  the  high  semblance  values  are  scatter¬ 
ing  of  higher  mode  surface  waves  to  P  waves  off  topographical  irregularities,  or  ran¬ 
dom  scattering  due  to  heterogeneity  in  the  area  of  the  test  site. 
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Introduction 

P  coda,  the  complex  waveform  that  follows  the  direct  P  arrival,  has  been  studied 
from  many  different  perspectives.  Coda  duration  has  been  used  for  many  years  to 
determine  the  size  of  local  earthquakes  (e.g.  Aid,  1969).  P  coda  amplitudes  have  also 
been  used  for  teleseismic  source  size  estimation,  because  it  is  generally  believed  that 
coda  consists  of  scattered  waves,  and  therefore  averages  out  the  effects  of  lateral 
heterogeneity  on  wave  amplitudes  by  sampling  a  large  volume  (Aki,  1982).  The 
scattering  models  have  led  to  widespread  use  of  coda  decay  as  a  method  for  determin¬ 
ing  attenuation  (e.g.  Aki  and  Chouet,  1975).  Several  workers  have  applied  random 
media  theory  (Chernov,  1960)  to  coda  in  order  to  define  the  strength  and  scale  length 
of  inhomogeneities  under  seismic  arrays  (Aki,  1973;  Capon,  1974;  Berteussen  et  al., 
1975).  Herraiz  and  Espinosa  (1987)  present  a  comprehensive  review  of  coda  research, 
and  their  terminology  will  be  followed  in  this  paper. 

While  there  has  been  substantial  success  in  explaining  various  properties  of  coda 
using  statistically  based  methods,  there  have  been  only  a  few  attempts  to  treat  coda  in 
a  deterministic  framework  by  identifying  actual  scattering  structures.  Key  (1967)  was 
able  to  distinguish  discrete  arrivals  in  the  P  coda  recorded  at  the  Eskdalemuir  array 
and  attributed  them  to  scattering  from  specific  topographic  features  near  the  array. 
Lay  (1987)  stacked  P  coda  from  underground  explosions  at  the  Pahute  Mesa  test  site 
in  Nevada  in  an  attempt  to  identify  contributions  from  isotropic  scatterers  in  the  source 
region.  However,  most  potential  scattering  locations  that  produced  any  waveform 
coherence  also  produced  enhanced  coherence  in  simulations  using  random  noise,  lead- 
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ing  to  some  ambiguity  in  the  final  result  In  this  study,  we  apply  the  basic  methodol¬ 
ogy  introduced  by  Lay  (1987)  to  explosions  at  the  nearby  Yucca  Flat  test  site  (Figure 
1),  in  an  attempt  to  identify  near-source  isotropic  scatterers. 

There  is  some  reason  to  anticipate  that  the  Yucca  Flat  source  array  might  be  a 
more  favorable  site  than  the  Pahute  Mesa  array  for  seeking  scattering  origins.  The 
Yucca  Flat  source  array  is  significantly  denser,  and  there  is  less  variation  in  explosion 
size.  Also,  a  study  by  Lynnes  and  Lay  (1988)  demonstrated  a  strong  correlation  of 
coda  excitation  with  position  in  the  Yucca  Flat  source  array,  indicating  significant 
near-source  contributions  to  the  P  coda. 

Methodology 

The  methodology  of  the  inversion  for  isotropic  scatterers  is  summarized  briefly 
below.  A  more  complete  development  is  given  in  Lay  (1987).  The  technique  essen¬ 
tially  consists  of:  assuming  an  x-y  location  where  primary  waves  traveling  at  horizon¬ 
tal  phase  velocity  vs  are  scattered  isotropically  to  teleseismic  P  waves;  computing  the 
appropriate  time  lags  for  each  source-station  pair,  stacking  the  seismograms  using 
these  time  lags;  and  computing  the  semblance  for  a  specified  time  gate.  A  high  sem¬ 
blance  value  for  a  given  x-y  location  and  assumed  moveout  velocity  (vs)  indicates  that 
a  coherent  arrival  can  be  identified,  suggesting  the  presence  of  a  structure  that  can  be 
approximated  as  an  isotropic  point  scatterer.  If  forward  or  back  scattering  dominate 
rather  than  isotropic  scattering,  the  coherence  will  diminish  but  should  not  be  elim¬ 
inated  so  long  as  many  seismograms  are  used. 
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The  scattered  arrival  time  lags  (Ty),  relative  to  the  direct  P  arrival,  are  given  by: 

xij  =  [(V*i)2  +  (yi“>i)2]V4K1  -  Pj  cos(0(H>y)3  (1) 

where  xit  yt  are  the  coordinates  of  the  source,  xr  ys  are  the  coordinates  of  the  hypothet¬ 
ical  scatterer,  pj  is  the  horizontal  ray  parameter  for  the  secondary  teleseismic  P  waves 
traveling  to  the  fh  station,  0,-  is  the  azimuth  from  the  source  to  the  hypothetical 
scatterer,  and  is  the  azimuth  from  the  scatterer  to  the  fh  station  (Figure  2).  The 
seismograms  are  then  stacked  after  applying  the  appropriate  time  lags,  and  the  sem¬ 
blance  (Neidell  and  Taner,  1971)  of  the  stack  is  computed  in  order  to  assess  whether 
any  waveform  coherence  can  be  associated  with  the  given  x—y  location  and  scattering 
velocity.  The  semblance  for  Ne  events  and  Nr  stations  is: 


*,){'  +  xij(.xsys,vs)) 


t  t=l  >=1 

where  Tg  is  the  time  gate  over  which  the  semblance  is  computed,  M  is  the  total 
number  of  seismograms,  s^t)  is  the  seismogram  from  the  ilh  source  at  the  fh  receiver, 
and  8,y  is  1  if  there  is  an  observation  for  the  ith  source  at  the  fh  receiver  and  0  if  there 
is  not.  The  semblance  can  be  viewed  as  the  power  of  the  stack  normalized  by  the 
total  power  in  the  seismograms,  and  it  varies  between  0  and  1,  unity  being  perfect 
waveform  coherence.  Semblance  is  also  closely  related  to  the  energy-normalized 
cross-correlation  (Neidell  and  Taner,  1971).  The  time  gate  ( Tg )  used  in  the  semblance 
computations  in  this  study  is  2  s,  corresponding  to  the  characteristic  pulse  duration  in 
the  short-period  recordings.  Tests  on  synthetic  data  using  a  longer  time  gate  of  3  s  did 
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not  significantly  change  the  semblances. 

The  F-test  can  be  applied  to  a  given  semblance  value  to  determine  whether  it  is 
statistically  significant.  The  F  value  for  a  non-central  distribution  (F')  is  related  to  the 
semblance  by  the  expression  (Douze  and  Laster,  1979): 

m/w* » o) 

N1  and  N2  arc  the  degrees  of  freedom,  and  k  is  the  non-centrality  parameter,  defined  as 
(Blandford,  1974): 


N{  =  2 BTg 

N2  =  IVj(M-l)  (4) 

Ps 

^ N 

where  B  is  the  bandwidth  of  the  signal,  P$  is  the  mean  power  of  the  signal  and  PN  is 
the  mean  power  of  the  noise  in  the  frequency  band  under  consideration.  The  F  statis¬ 
tic  for  a  non-central  distribution  (F)  can  be  approximated  by  the  F  statistic  for  a  cen¬ 
tral  distribution  (Zelen  and  Severo,  1965): 


F(F'  |  Nlt  N2,  k)  =  F(F  |  Nx*,  N?) 
Nx  p' 


F  =  F- 


Nx  +  k  1+F S/PN 


(5) 


(Nx+k)2 

“  (Nx+2k)  ~Nl 


1  + 


(W/)2 


(1  +  2F 5/P N) 


The  statistical  significance  can  then  be  assessed  from  standard  tables  using  the  F  value 
with  Nx*  and  N2  degrees  of  freedom.  The  validity  of  the  statistical  measure  requires 
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assumptions  about  the  noise  properties,  which  are  not  generally  satisfied  by  the  raw 
data.  However,  weighting  functions,  such  as  those  described  in  the  data  analysis  sec¬ 
tion  below,  can  be  used  to  suppress  correlated  noise,  allowing  the  use  of  the  F-test. 

Data  Analysis. 

The  P  waves  that  we  analyze  here  were  recorded  on  short-period  instruments  at 
65  teleseismic  (25°-95°)  WWSSN  and  Canadian  Seismic  Network  (CSN)  stations. 
Records  from  32  underground  explosions  at  the  Yucca  Flat  test  site  were  used,  yield¬ 
ing  a  total  of  1050  recordings.  The  first  25  s  of  the  P  wavetrain  were  digitized  on  all 
of  the  seismograms.  In  order  to  obtain  a  uniform  database,  the  CSN  records  were 
multiplied  in  the  frequency  domain  by  the  ratio  of  the  standard  WWSSN  instrument 
response  to  the  individual  CSN  instrument  responses,  thus  equalizing  the  CSN  seismo¬ 
grams  to  WWSSN- type  recordings.  A  zero-phase  three-pole  Butterworth  bandpass 
filter  from  0.2  to  2.0  Hz  was  applied  to  minimize  digitization,  instrument,  and  earth 
noise.  The  seismograms  were  also  normalized  to  the  peak  amplitude  in  order  to 
suppress  the  amplitude  effects  of  variable  explosion  yield,  attenuation  and 
focusing/defocusing.  In  order  to  obtain  more  uniform  azimuthal  coverage,  the  seismo¬ 
grams  within  7.5°  azimuth  windows  were  averaged,  reducing  the  number  of  indepen¬ 
dent  observations  from  1050  to  816. 

The  semblance  procedure  relies  on  correlating  similar  features  occurring  at 
different  times  in  the  full  set  of  seismograms  due  to  the  directivity  of  the  scatterei. 
Consequently,  it  is  helpful  to  remove  features  that  occur  at  the  same  time  in  a  large 
subset  of  seismograms.  The  most  obvious  such  feature  is  the  direct  arrival  itself, 
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which  normally  has  distinctive  characteristics  due  to  both  source  and  receiver  effects. 
Other  coherent  arrivals  might  be  produced  by  slapdown  or  aftershocks,  which  would 
be  common  to  all  seismograms  for  a  given  source,  and  scattered  arrivals  from  the 
vicinity  of  a  given  station,  which  would  be  similar  for  all  seismograms  at  that  receiver. 
In  an  attempt  to  minimize  the  effects  of  these  features,  time-dependent  weighting 
filters  were  applied  to  the  seismograms  (s^t))  before  stacking: 

s'ijit)  =*$)[  l-wvr/r)]  [  l-wy/r)]  (6) 

where  we  and  wr  are  constants  between  0  and  1 ,  e,{r)  is  the  normalized  envelope  of  the 
stacked  traces  for  the  ***  event,  and  rj(t)  is  the  normalized  envelope  of  the  stacked 
traces  for  the  fh  station.  These  filters  downweight  receiver  effects  and  the  direct 
arrival,  which  would  otherwise  dominate  the  semblance  computations  and  result  in 
high  semblances  centered  on  the  sources.  A  wide  range  of  values  and  combinations  of 
the  weighting  constants  we  and  wr  was  explored,  with  a  value  of  0.75  for  both  con¬ 
stants  producing  good  results  in  simulations  using  synthetic  data  sets.  This  choice  of 
we  and  wr  is  kept  fixed  below. 

Tests  Using  Synthetic  Data 

In  order  to  assess  the  intrinsic  ability  of  our  particular  source-station  geometry  to 
resolve  the  locations  and  scattering  velocities  of  isotropic  point  scatterers,  we  con¬ 
ducted  several  synthetic  experiments.  Synthetic  seismograms  were  computed  compris¬ 
ing  a  direct  arrival  and  a  P  wave  scattered  at  a  specified  moveout  velocity  and  loca¬ 
tion,  with  a  prescribed  amplitude  relative  to  the  direct  arrival.  The  seismograms  were 
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formed  using  a  modified  Haskell  source  (Helmberger  and  Hadley,  1981),  a  short-period 
WWSSN  instrument  response  and  a  Futterman  Q  operator  with  t*  of  0.4  s.  A  seismo¬ 
gram  was  constructed  for  each  source-station  pair  contributing  to  the  actual  data  set 
The  semblance  was  then  calculated  at  each  point  on  a  square  19x19  grid  centered  on 
the  Yucca  Flat  test  site,  with  a  grid  spacing  of  2  km.  The  x  coordinate  increases  to  the 
east  and  the  y  coordinate  increases  to  the  north.  The  results  of  a  simulation  for  a 
point  scatterer  with  a  horizontal  phase  velocity  ( vs )  of  2.5  km/s  are  shown  in  Figure  3. 
The  figure  on  the  left  (Figure  3a)  shows  the  grid  points  at  which  semblances  were 
computed,  with  the  closed  circle  indicating  the  position  of  the  scatterer  producing  the 
"coda"  arrival  in  the  synthetic  seismograms.  A  perspective  plot  of  the  semblance  at 
each  grid  point  computed  for  this  seismogram  suite,  assuming  a  scattering  velocity  of 
2.5  km/s,  is  shown  in  Figure  3b.  The  highest  semblance  (1.0)  is  found  at  the  correct 
location  for  the  scatterer,  although  several  neighboring  locations  also  produce  high 
semblance  values  due  to  spurious  coherence,  resulting  in  limited  spatial  sensitivity,  or 
"resolution".  (Note  that  this  is  not  a  standard  mathematical  inversion  in  that  the  sem¬ 
blance  at  a  given  point  is  independent  of  the  semblance  at  other  points,  and  model 
parameters  are  not  obtained  directly.  Hence,  the  resolution  is  not  strictly  defined,  and 
behaves  differently  than  for  a  linear  inversion.) 

We  also  attempted  to  assess  the  effect  of  noise  on  the  semblance  calculations  by 
adding  arrivals  with  random  delay  times  and  amplitudes  to  the  synthetic  seismograms 
for  a  single  point  scatterer.  The  effect  of  adding  50  random  arrivals  to  each  waveform 
in  the  previous  synthetic  seismogram  set  is  shown  in  Figure  3c.  The  random  arrivals 
have  uniform  distributions  of  amplitude  relative  to  the  direct  arrival  (with  ratios 
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between  -0.5  and  0.5),  and  of  arrival  times  within  the  window  of  the  seismogram. 
The  peak  semblance  is  significantly  reduced,  as  might  be  expected  due  to  the  addition 
of  noise.  However,  the  noise  reduces  the  semblance  values  of  the  surrounding  grid 
points  even  more,  resulting  in  increased  spatial  sensitivity  for  the  synthetic  data  set. 

At  first  sight,  the  apparent  enhancement  of  "resolution"  with  the  addition  of  noise 
is  somewhat  paradoxical.  The  explanation  lies  in  the  nature  of  the  semblance  statistic. 
Since  semblance  is  formed  by  summing  lagged  traces  rather  than  multiplying  them, 
small  misalignments  can  still  produce  high  semblance  values  if  the  waveforms  are 
sufficiently  similar.  When  random  arrivals  arc  added  to  the  waveforms,  they  degrade 
the  semblance  at  all  of  the  locations.  However,  the  large  number  of  waveforms  still 
produces  a  significant  semblance  at  the  precise  location  for  which  the  waveforms  are 
correctly  aligned. 

Another  apparent  paradox  in  the  high  resolution  is  that  the  time  differences  asso¬ 
ciated  with  nearby  grid  points  are  quite  small,  about  0. 1-0.2  s.  Thus  the  resolution  is 
better  than  might  be  expected  from  the  high  frequency  comer  of  the  filter  (2  Hz).  This 
is  due  to  the  influence  of  phase  information  in  the  semblance  technique,  which  can  be 
recovered  even  from  the  lower  frequencies. 

The  ability  of  the  technique  to  resolve  the  horizontal  phase  velocity  of  the  waves 
scattered  by  a  point  scatterer  was  tested  by  computing  the  semblances  for  several 
different  assumed  moveout  ve’ocities,  using  the  same  suite  of  synthetics  as  the  above 
simulation  of  a  point  scatterer  plus  random  noise  (Figure  4).  The  highest  peak  sem¬ 
blance  value  occurs  for  the  semblances  calculated  for  a  velocity  of  2.5  km/s,  the 
correct  velocity  (Figure  4b).  The  low  semblances  for  the  other  scattering  velocities 
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(Figure  4a,  c)  indicate  that  the  velocity  can  be  resolved  well  for  a  single  point  scatterer 
in  this  source- station  geometry.  This  is  important  for  identifying  what  type  of  scatter¬ 
ing  is  occurring. 

Since  the  variation  of  delay  times  of  scattered  arrivals  depends  strongly  on  the 
source-station  geometry,  the  peak  semblance  and  spatial  resolution  may  vary  in 
different  parts  of  the  source  region.  To  characterize  these  variations,  we  computed  81 
suites  of  synthetic  seismograms,  assuming  single  point  scatterers  at  x,y  (E,  N)  =  0  km, 
±4  km,  ±8  km,  ±12  km,  and  ±16  km  (Figure  5),  scattering  waves  with  a  scattering 
velocity  (v,)  of  2.5  km/s.  The  point  scatterers  radiated  waves  with  half  the  amplitude 
of  the  direct  arrival,  and  50  random  noise  arrivals  were  also  added  to  each  seismo¬ 
gram.  The  semblances  for  the  full  19x19  grid  (x=-18  — »  +18  km,  y=-18  — >  +18  km) 
were  then  calculated  for  each  set  of  seismograms  at  the  correct  stacking  velocity  (2.5 
km/s).  The  semblance  surfaces  for  four  of  the  simulations  are  shown  in  Figure  6.  In 
each  case,  the  actual  location  of  the  scatterer  can  be  identified  as  the  point  with  the 
highest  semblance.  Interestingly,  the  peak  semblance  for  a  point  scatterer  in  the  mid¬ 
dle  of  the  array  (Figure  6d)  is  smaller  than  for  scatterers  in  other  locations.  This  is 
due  to  contamination  from  the  direct  arrivals  and  the  effects  of  the  weighting  func¬ 
tions.  Another  interesting  feature  is  the  anisotropic  nature  of  the  resolution  for  points 
outside  the  array  (Figure  6a-c),  with  the  semblances  dropping  off  from  the  peak  sem¬ 
blance  less  rapidly  in  the  NE-SW  direction  than  in  the  NW-SE  direction. 

The  peak  semblances  for  all  81  simulations  are  shown  in  Figure  7.  As  mentioned 
above,  point  scatterers  near  the  middle  of  the  source  array  actually  produce  relatively 
low  semblance  values.  Excluding  these  locations,  the  peak  semblances  are  roughly 
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equal  for  different  positions  in  the  source  region. 

We  chose  to  characterize  the  variable  spatial  resolution  by  a  two-dimensional 
"spread  function".  The  spread  function  was  computed  by  comparing  the  semblance  at 
an  actual  scatterer  location  (5max)  with  the  semblances  calculated  at  each  of  the  eight 
adjacent  points.  For  each  azimuth  (0)  from  the  actual  location  to  the  adjacent  point, 
the  spread  (cr(0))  is  defined  in  terms  of  5^  and  the  semblance  at  the  adjacent  point  in 
that  direction,  5,  as: 


0(0)  =  1  - 


AS 


'max 


AS  = 


AS  =  S^-S  for  0=0°.  90°.  180°,  and  270°, 

(‘-'max- *S") 


(7) 


V2 


for  0=45°,  135°,  225°,  and  315°. 


The  spread  in  a  given  direction  is  equal  to  1  for  negligible  spatial  resolution  over  the 
given  grid  spacing  (SmaT=S)  and  is  equal  to  zero  for  perfect  resolution,  (5=0).  Thus, 
the  spread  function  is  essentially  reciprocal  to  the  resolution.  The  results  are  presented 
as  octagons,  with  the  distance  from  the  center  to  each  vertex  proportional  to  the  spread 
function  in  that  direction.  A  large  octagon  represents  a  large  spread  function  in  all 
directions,  implying  poor  spatial  resolution  of  a  scatterer  at  that  point,  while  a  small 
octagon  indicates  a  small  spread  function  and  good  resolution.  An  equidimensional 
octagon  would  indicate  equal  resolution  in  each  direction.  However,  most  of  the 
spread  functions  show  a  pronounced  elongation  in  the  NE-SW  direction,  indicating 
poorer  resolution  in  this  direction  t!,an  in  other  directions.  This  is  probably  due  to  the 
roughly  linear  NW-SE  trend  of  the  source  array.  The  best  resolution  (smallest  spread 
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function)  appears  to  be  immediately  to  the  ENE  and  WSW  of  the  array.  These  charac¬ 
teristics  can  be  compared  with  the  analogous  variations  for  the  Pahute  Mesa  array 
geometry  examined  by  Lay  (1987)  to  gain  an  appreciation  for  the  effects  of  source  dis¬ 
tribution. 


Inversion  of  Real  Data 

The  actual  data  set  was  inverted  for  isotropic  point  scatterers  using  the  same  grid 
as  in  the  synthetic  cases  (Figures  4-7),  for  moveout  velocities  of  1.5  km/s,  2.5  km/s 
and  3.5  km/s  (Figure  8).  All  of  the  semblance  surfaces  have  relatively  high  values  in 
the  center.  The  two  higher  moveout  velocities  also  appear  to  have  a  high  semblance 
area  about  10  km  to  the  ENE  of  the  Yucca  Flat  source  array,  and  a  ridge  of  moderate 
values  10  km  WNW  of  the  array.  Both  of  these  areas  are  regions  of  fairly  good  reso¬ 
lution,  according  to  the  simulations  (Figure  7).  However,  the  high  semblance  area  in 
the  data  does  seem  to  be  somewhat  elongated  in  the  direction  orthogonal  to  the  trend 
of  the  source  array  (i.e.,  NE-SW),  as  predicted  by  the  resolution  experiment  (Figure  7). 
Semblances  were  also  computed  for  moveout  velocities  of  1.0  km/s,  2.0  km/s,  3.0 
km/s,  and  4.0  km/s.  The  high  semblance  feature  to  the  ENE  appears  for  all  velocities 
higher  than  2.0  km/s.  All  stacking  velocities  produce  high  semblances  in  the  central 
region. 

In  order  to  test  the  stability  of  the  results,  the  semblance  was  computed  for  three 
subsets  of  the  data:  the  16  most  northerly  events,  the  16  most  southerly  events,  and  the 
16  largest  events.  In  each  case,  the  overall  semblance  values  for  all  of  the  moveout 
velocities  were  higher  by  a  factor  of  2-3  for  the  subsets  than  for  the  full  set  of  data. 
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Qualitatively,  the  semblance  surfaces  for  the  subsets  were  quite  similar  to  those  for  the 
full  set,  except  that  the  semblances  for  the  subsets  were  disproportionately  high  in  the 
center  of  the  region.  This  is  probably  due  to  spurious  waveform  coherence  associated 
with  the  direct  arrivals. 

Unlike  the  simulations,  the  peak  semblances  for  the  full  set  of  real  data  are 
extremely  low,  on  the  order  of  0.02-0.03  (Figure  8).  However,  these  values  are  actu¬ 
ally  somewhat  higher  than  those  obtained  for  Pahute  Mesa  by  Lay  (1987)  when  station 
and  event  weighting  were  applied.  In  fact,  a  semblance  of  0.02  yields  an  F'  value  of 
16.6,  with  degrees  of  freedom  A^=7.2  and  N2=5 868.  Of  course,  the  statistical 
significance  depends  on  the  signal  to  noise  ratio,  which  is  unknown.  If  we  assume  a 
signal/noise  ratio  of  0,  then  F=F',  and  a  semblance  of  0.02  is  significant  at  the  99% 
level.  On  the  other  hand,  if  we  assume  that  the  signal/noise  power  ratio  is  equal  to  2, 
F=F'/3=5.5,  and  N{*  =  1.8Ar1  =  13.0,  which  is  still  significant  at  the  99%  level. 

In  any  case,  even  the  highest  semblances  in  the  real  data  are  low  enough  to  cause 
concern  that  they  might  be  due  simply  to  incomplete  suppression  of  noise  caused  by 
the  source-station  geometry  To  rest  this  possibility,  semblances  were  computed  for  a 
suite  of  synthetic  seismograms  comprising  the  direct  arrival  and  50  random  arrivals, 
varying  in  amplitude  between  -0.5  and  0.5  relative  to  the  direct  arrival.  The  synthetic 
seismograms  were  computed  for  the  same  source-station  pairs  as  the  data.  The  sem¬ 
blances  were  computed  assuming  a  stacking  velocity  of  2.5  km/s  (Figure  9a).  The 
moveout  velocity  of  2.5  km/s  was  chosen  because  it  appears  to  produce  more  stable 
semblance  surfaces  for  the  data  than  the  higher  moveout  velocities.  A  velocity  of  2.5 
km/s  is  probably  too  fast  for  short-period  fundamental  mode  Rayleigh  waves,  since  the 
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S  wave  velocities  of  the  basin  sediments  are  about  1  km/s  and  the  Paleozoic  basement 
is  about  2.6  km/s  (Ferguson,  1982).  However,  this  is  a  reasonable  velocity  for  higher 
mode  surface  waves. 

The  peak  semblance  in  the  random  simulation  is  surprisingly  high  (0.11),  suggest¬ 
ing  that  many  of  the  higher  semblances  in  the  data  stacks  may  be  artifacts.  However, 
the  high  semblances  are  mostly  confined  to  the  center  of  the  grid,  in  proximity  to  the 
source  array.  This  is  likely  due  to  the  residual  coherence  of  the  direct  arrivals,  which 
is  not  completely  removed  by  the  station  and  event  weighting  filters.  The  lack  of 
appreciable  semblance  in  the  random  simulations  at  outlying  locations  suggests  that 
random  noise  can  not  explain  high  semblance  in  the  data  (Figure  9b)  away  from  the 
central  region  of  the  source  array.  The  random  noise  simulations  produce  simpler 
semblance  surfaces  than  do  similar  simulations  for  Pahute  Mesa  (Lay,  1987),  probably 
due  to  the  closer  spacing  of  the  Yucca  Flat  events. 

Since  the  highest  semblance  values  occur  in  a  region  where  the  spatial  resolution 
is  greatest,  it  is  also  possible  that  the  distribution  of  scatterers  is  random,  but  high 
semblances  are  found  only  for  points  where  the  spatial  resolution  is  favorable.  In 
order  to  test  this  hypothesis,  we  computed  synthetic  seismograms  for  a  set  of  70 
scatterers  with  random  amplitudes  up  to  0.2  the  amplitude  of  direct  arrival,  stacking 
velocity  of  2.5  km/s,  and  positions  within  10  km  of  the  center  of  the  grid.  No  random 
noise  arrivals  were  added  to  the  seismograms.  The  locations  and  amplitudes  of  these 
scatterers  are  shown  in  Figure  9c.  The  semblance  was  then  computed  for  each  grid 
point  assuming  a  scattering  velocity  of  2.5  km/s  (Figure  9d).  The  semblances  are 
again  quite  high  (up  to  0.24).  In  contrast  to  the  random  noise  case,  however,  the 
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random  scatterer  simulation  produces  relatively  low  semblance  values  in  the  central 
region,  but  higher  semblances  in  the  area  to  the  northeast  where  the  high  semblances 
were  located  in  the  data  inversions.  Thus,  while  the  technique  finds  high  semblance  in 
some  areas  where  scatterers  are  located,  it  fails  in  this  case  to  locate  other  scatterers, 
apparently  due  to  the  sensitivity  characteristics  imposed  by  the  array  geometry. 

Discussion 

Several  studies  have  indicated  that  topography  is  a  major  contributor  to  P  wave 
coda.  Using  a  statistical  approach,  Greenfield  (1971)  was  able  to  model  the  energy 
flux  of  coda  from  Novaya  Zemlya  explosions  by  hypothesizing  Rayleigh— >P  scattering 
by  the  island  coastline,  river  valleys  and  other  topography.  Key  (1967)  used  standard 
array  analysis  to  determine  the  location  and  phase  velocity  of  coherent  arrivals  in  the  P 
coda  at  the  Eskdalemuir  array.  The  locations  were  generally  associated  with  topo¬ 
graphic  highs.  Key  (1967)  found  one  particularly  strong  arrival  to  be  associated  with  a 
river  valley,  with  a  phase  velocity  of  2.5  kirVs,  suggesting  P— »RayIeigh  scattering  near 
the  receiver.  Synthetic  calculations  also  indicate  a  scattering  role  for  topographical 
highs.  McLaughlin  and  Jih  (1988)  computed  2-D  finite  difference  synthetics  for  Ray¬ 
leigh  waves  incident  on  a  mesa.  The  resulting  wavefield  includes  both  forward  and 
back  scattered  P  waves.  In  order  to  see  if  topography  can  explain  the  high  semblance 
regions  located  in  this  study,  the  semblances  computed  from  the  data  with  a  stacking 
velocity  of  2.5  km/s  are  superimposed  on  a  topographic  map  of  the  region  in  Figure 
10.  Indeed,  the  isolated  high-semblance  area  to  the  east  is  associated  with  a  similarly 
isolated  topographic  high  in  the  Halfpint  Range.  The  ridge  to  the  west  also  correlates 
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with  moderately  high  semblance  values.  Thus,  waves  scattered  from  topographical 
irregularities  seem  to  be  a  plausible  explanation  for  some  of  the  early  P  coda  for 
Yucca  Flat  explosions.  However,  the  semblances  in  the  northwest  area  of  the  grid, 
where  the  topography  is  the  highest,  are  rather  low.  If  the  topography  does  produce 
the  majority  of  the  coda,  the  low  semblances  might  be  due  to  the  shape  of  the  topo¬ 
graphic  high.  On  the  other  hand,  the  random  scatterer  simulation  indicates  that  an 
equally  valid  explanation  of  the  semblances  derived  from  the  data  would  be  a  random 
distribution  of  scatterers  in  the  central  region  of  the  test  site,  which  is  quite  flat 

In  summary,  the  semblance  technique  has  yielded  somewhat  ambiguous  results  in 
terms  of  identifying  deterministic  scatterers.  Indeed,  the  simulations  seem  to  indicate 
that  random  scatterers  explain  the  data  equally  well.  The  main  problem  seems  to  be 
that  the  source  array  geometry  exercises  strong  control  on  the  sensitivity  of  the  tech¬ 
nique.  This  has  been  the  case  both  for  the  two-dimensional  Pahute  Mesa  source  array 
(Lay,  1987)  as  well  as  the  linear  Yucca  Flat  array. 

However,  an  additional  factor  may  be  the  limitations  of  the  data  set,  which  may 
obscure  more  significant  features  that  would  yield  high  semblance  in  areas  of  the  test 
site  that  are  not  favored  by  the  array  geometry.  For  instance,  the  use  of  analog,  rather 
than  digital,  data  places  a  severe  constraint  on  the  frequency  range  that  can  be  reliably 
used.  Digital  data  bases  of  substantial  size  are  now  accumulating  for  more  recent 
events.  Another  factor  that  probably  decreases  the  semblance  of  the  waveforms  is  the 
influence  of  lateral  heterogeneity  on  the  travel  time  from  the  source  to  the  scatterer. 
Accounting  for  this  effect  could  result  in  improved  alignment  of  the  waveforms.  Also, 
isotropic  scattering  was  assumed,  a  condition  that  could  be  relaxed  to  allow  for 
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radiation  patterns.  For  NTS,  this  assumption  is  not  particularly  unrealistic:  due  to  the 
low  velocities  near  the  surface,  the  teleseismic  P  waves  leave  the  source  area  in  a  very 
naiTOw  cone  of  rays.  Furthermore,  the  synthetic  results  of  McLaughlin  and  Jih  (1988) 
indicate  that  Rayleigh— scattering  from  topography  is  roughly  isotropic  for 
wavelengths  that  are  significantly  greater  than  the  height  of  the  irregularity.  One 
modification  of  the  technique  that  might  produce  some  improvement  in  the  results 
would  be  an  inversion  of  the  coda  envelopes,  rather  than  the  seismograms  themselves. 
This  would  mitigate  the  effects  of  phase  incoherence  and  receiver  and  source  charac¬ 
teristics  that  influence  the  dominant  frequencies  in  the  seismograms.  However,  the  use 
of  envelopes  would  dramatically  decrease  the  bandwidth;  also,  the  statistics  of  cross- 
correlation  of  such  one-sided  functions  are  not  so  well-developed  as  for  oscillatory 
waveforms. 


Conclusion 

P  wave  seismograms  from  32  Yucca  Flat  explosions  were  inverted  for  isotropic 
point  scatterers  in  the  source  region  by  computing  the  semblance  for  hypothetical 
scatterers  at  a  series  of  locations.  Synthetic  experiments  show  the  ability  of  the  sem¬ 
blance  analysis  to  identify  strong  point  scatterers,  with  the  attainable  resolution  varying 
with  position  and  direction  in  the  source  region.  The  resolution  is  generally  best 
within  a  few  kilometers  of  the  source  array,  and  is  better  in  the  NW-SE  direction  than 
in  the  NE-SW  direction  due  to  the  trend  of  the  source  array. 

The  semblances  of  the  real  data,  while  indicating  much  lower  waveform  coher¬ 
ence  than  for  the  simulations,  are  relatively  high  about  10  km  ENE  of  the  source  array, 
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and  moderately  high  about  10  km  WNW  of  the  array,  as  well  as  in  the  immediate 
vicinity  of  the  array,  for  a  moveout  velocity  of  25  km/s.  Analysis  of  signals  with  ran¬ 
dom  noise  yields  high  semblances  in  the  center  as  well,  suggesting  that  semblances  in 
this  area  may  be  artifacts  of  the  source-station  geometry  and  contamination  by  the 
direct  arrivals.  However,  random  noise  does  not  appear  to  explain  the  high  semblance 
values  to  the  ENE  and  WNW.  The  high  semblance  area  to  the  ENE  is  associated  with 
an  isolated  topographic  high,  suggesting  that  topographic  irregularities  may  play  an 
important  role  in  coda  generation.  However,  a  simulation  of  random  scatterers  pro¬ 
duced  high  semblances  in  the  same  area  ENE  of  the  array  as  the  data,  indicating  that  a 
random  distribution  of  scatterers  in  the  vicinity  of  the  test  site  might  also  explain  the 
results. 
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Tigure  3  (a)  Locations  for  which  semblance  was  calculated  The  diamonds  show' 

the  locations  of  the  Yucca  Flat  underground  explosions  used  in  this  studv 
The  closed  circle  is  the  location  of  the  point  scatterer  contributing  to  the 
suite  of  synthetic  seismograms,  (b)  Perspective  plot  of  semblances  for  the 
corresponding  suite  of  synthetic  seismograms.  The  peak  semblance  value  is 
shown  in  the  upper  left,  let  Perspective  plot  of  semblances  for  a  suite  of 
seismograms  comprising  .ltn'als  from  the  scatterer  in  (a),  plus  50  random 
arrivals  in  each  seismogram 
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Figure  4.  Semblances  computed  for  the  suite  of  synthetic  seismograms  for  the 
scatterer  in  Figure  3a,  plus  random  noise.  Three  different  stacking  velocities 
are  shown,  with  the  actual  one  in  the  middle.  The  true  position  of  the 
scatterer  is  shown  by  the  circle,  which  is  closed  for  the  correct  velocity  and 
open  for  the  incorrect  velocities. 
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Figure  5.  Locations  of  the  point  scatterers  used  in  81  single  point  scatterer 
simulations,  together  with  the  grid  of  locations  used  in  the  semblance 
analysis  of  the  synthetic  data  sets. 
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Figure  6.  Semblances  computed  for  point  scatteres  at  four  of  the  locations  in  Figure 
5.  The  numbers  in  the  upper  left  are  the  peak  semblances  tor  each 
simulation,  and  the  numbers  in  the  upper  right  are  the  x-y  coordinates  of 
the  point  scatterer  (located  at  the  closed  circle). 
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Figure  7.  (a)  Semblance  values  (Smax)  at  the  actual  scatterer  locations  for  the  81 

single  point  scatterer  simulations  (see  Figure  5).  (b)  Octagons  showing 

spread  functions  in  eight  directions  (N,  S,  E,  W,  NE,  NW,  SE  and  SW)  for 
the  same  81  simulations,  indicated  by  the  radius  to  each  vertex  of  the 
octagon.  The  spread  function  varies  between  0  (perfect  spatial  resolution 
near  the  peak  semblance)  and  1  (poor  spatial  resolution).  An 
equidimensional  octagon  indicates  equal  spatial  resolution  in  all  directions, 
and  the  best  overall  resolution  is  for  the  smaller  symbols,  oee  text  for  a 
more  comprehensive  explanation. 
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Figure  8.  Semblances  calculated  for  the  Yucca  Flat  data  set,  assuming  three  different 
scattering  velocities.  The  peak  semblance  for  each  surface  is  shown  in  the 
upper  left. 
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max.  amplitude:  0.2  max.  semblance=0.24 

Figure  9.  (a)  Map  view  of  semblances  calculated  for  a  random  simulation,  (b) 
Map  view  of  semblances  from  data,  (c)  Map  view  of  randomly  distributed 

scatterers  used  in  random  scatterer  simulation.  The  size  of  the  symbol  indicates 
the  amplitude  of  the  scattered  arrivals  (relative  to  the  direct  arrival,  (d)  Map 

view  of  random  scatterer  simulation.  Note  the  similarity  to  the  map  view  of 
the  data  semblances.  The  peak  value  for  each  surface  is  shown  below  each 
graph,  and  a  stacking  velocity  of  2.5  km/s  was  used  for  all. 


Section  4 

Effects  of  Lateral  Velocity  Heterogeneity  Under  the  Nevada  Test  Site  on 
Short  Period  P  Wave  Amplitudes  and  Travel  Times 

(Accepted  for  publication  in  Pure  and  Applied  Geophysics;  authors  Christopher 
S.  Lynnes  and  Thorne  Lay) 
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ABSTRACT 

Short-period  teleselsmic  P  waves  from  the  Nevada  Test  Site  (NTS)  show  systematic 
variations  in  amplitudes  and  travel  times,  with  low  amplitudes  corresponding  to  fas: 
travel  times,  suggesting  elastic  focussing-defocussing  effects.  Also,  the  patterns  for 
events  at  the  Pahute  Mesa  subsite  are  systematically  different  from  those  at  the  Yucca 
Flat  subsite,  indicating  the  presence  of  a  near-source  component  in  both  the  amplitude 
and  travel  time  variations.  This  component  is  isolated  by  removing  the  mean  station 
pattern  for  all  of  NTS  from  the  observations.  A  very-near-source  component  in  the 
Pahute  Mesa  observations  is  also  isolated  by  removing  subsite  station  means  from  the 
measurements.  These  anomalies  are  back-projected  through  laterally  homogeneous 
structure  to  form  thin  lens  models  at  various  depths.  Travel  time  delays  are  predicted 
from  the  amplitude  variations  using  the  equation  for  wavefront  curvature.  The  long- 
wavelength  components  of  the  predicted  and  observed  time  delays  correlate  well  at 
depths  of  25  km  for  the  very-near-source  component  under  Pahute  Mesa  and  160  km 
for  the  regional  component  under  NTS.  The  time  delay  surfaces  predicted  by  the 
amplitudes  at  these  depths  are  mapped  into  warped  velocity  discontinuities  suitable  for 
the  calculation  of  synthetic  seismograms  using  the  Kirchhoff- Helmholtz  integral 
formulation.  Both  the  intersite  (near-source)  and  intrasite  (very-near-source) 
differences  in  amplitudes  are  qualitatively  predicted  very  well,  although  the  range  of 
variation  is  somewhat  underpredicted.  This  deficiency  is  likely  due  to  the  destructive 
interference  of  anomalies  inherent  in  back-projection  to  a  single  layer. 
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INTRODUCTION 

Theoretically,  an  underground  explosion  should  radiate  compressional  waves 
isotropically.  However,  several  studies  of  short-period  P  waves  from  explosions  (e.  g.. 
Lay  et  al.,  1984;  Burger  et  al.,  1986)  have  indicated  that  amplitudes  can  have  coherent 
variations  with  azimuth  ranging  over  factors  of  3  to  10.  Such  variations  are  important  to 
understand  because  of  the  uncertainty  they  add  to  underground  test  size  estimation. 
Furthermore,  those  variations  can  be  extremely  useful  in  investigating  lateral  variations 
in  earth  structure,  particularly  when  examined  in  concert  with  travel  time  data.  Given 
sufficiently  accurate  earth  models,  it  may  be  possible  to  correct  the  amplitude  variations 
for  propagational  effects,  thereby  improving  the  accuracy  of  the  size  estimates. 

In  this  study,  we  examine  short-period  P  wave  amplitudes  in  conjunction  with  travel 
times  from  underground  explosions  at  the  Nevada  Test  Site  (NTS).  A  previous  analysis 
(Lynnes  and  Lay,  1987)  has  established  that  low  amplitudes  correlate  with  fast  travel 
times,  implying  that  variations  in  velocity,  rather  than  attenuation,  are  primarily 
responsible  for  the  amplitude  fluctuations.  These  amplitude  variations  are  used  in  this 
study  to  derive  laterally  varying  velocity  structures,  parameterized  as  warped  velocity 
discontinuities  that  are  as  compatible  as  possible  with  the  travel  time  observations. 

Underground  explosions  at  NTS  (Fig.  1)  provide  a  dense,  well-located  array  of 
sources  that  offers  an  excellent  site  for  such  a  study.  Also,  several  travel  time  studies 
have  already  found  evidence  for  significant  lateral  heterogeneity  in  the  crust  and  upper 
mantle  beneath  NTS  (Spence,  1974;  Minster  et  al.,  1981;  Taylor,  1983),  and  the  geology 
of  the  general  area  >s  replete  with  Cenozoic  calderas,  which  might  be  expected  to 
produce  or  be  related  to  lateral  velocity  variations. 

A  few  previous  studies  have  predicted  amplitude  variations  using  models  derived 
from  travel  time  data.  Thomson  and  Gubbins  (1982)  calculated  amplitudes  for  the 
Norwegian  Seismic  Array  (NORSAR)  using  a  smoothly  varying  structure  obtained  by  a 
linear  inversion  of  travel  times.  Although  they  found  statistically  significant  correlations 
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between  predicted  and  observed  amplitudes,  the  predicted  amplitude  range  was  much 
smaller  than  the  observed  range.  Cormier  (1987)  traced  rays  through  Taylor’s  (1983) 
model  of  the  lithosphere  under  NTS,  producing  amplitude  variations  for  a  Pahute  Mesa 
explosion  that  were  similar  to  the  observed  variations  reported  in  Lay  et  al.  (1984).  On 
the  other  hand.  Lay  et  al.  (1986)  traced  rays  through  the  structure  derived  by  Minster  et 
al.  (1981)  from  travel  times  for  the  same  area,  finding  almost  no  variation  in  amplitudes. 
The  two  travel  time  inversions  use  overlapping  data  sets  of  teleseismic  travel  times,  but 
differ  in  the  parameterization  of  the  model:  Taylor’s  model  extends  to  100  km  depth, 
and  includes  information  from  Pn  source  terms,  while  Minster  et  al.  invert  for 
perturbations  down  to  ISO  km  and  employ  an  additional  smoothing  algorithm  in  the 
inversion.  The  result  is  that  Taylor’s  (1983)  model  has  velocity  perturbations  that  are 
larger  than  in  Minster  et  al.'s  model  by  a  factor  of  three  and  fluctuate  more  rapidly. 

The  attendant  differences  in  lateral  velocity  gradients  between  the  two  models  have  a 
drastic  effect  on  the  predicted  amplitude  variations. 

Given  the  dependence  of  predicted  amplitudes  on  the  parameterization  of  the  travel 
time  models,  it  would  be  desirable  to  obtain  a  velocity  structure  from  the  amplitudes 
themselves.  However,  due  to  the  non-linear  relationship  of  amplitudes  to  velocity 
structure,  such  inversions  have  been  few.  Haddon  and  Husebye  (1978)  projected 
amplitude  variations  onto  surfaces  at  depths  between  150  and  200  km  below  the 
NORSAR  array,  and  used  an  equation  describing  wavefront  curvature  to  predict  travel 
time  variations  due  to  an  infinitely  thin  lens,  obtaining  a  high  correlation  with  observed 
travel  time  residuals.  Thomson  (1983)  inverted  amplitudes  recorded  at  NORSAR  for  a 
smoothly  varying  structure  using  a  raytracing  formalism,  obtaining  fairly  good  amplitude 
predictions  but  only  moderate  correlations  of  predicted  and  observed  travel  times.  Thus, 
amplitudes  seem  to  contain  information  that  is  partially  complementary  to  that  in  the 


travel  times. 
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Another  attempt  to  produce  the  amplitude  variations  at  NTS  is  the  forward  modeling 
approach  of  Scott  and  Helmberger  (1985),  who  assumed  an  up  warped  Mohorovicic 
discontinuity.  By  calculating  synthetic  seismograms  for  arbitrary  structures  using  the 
Kirchhoff-Helmholtz  integral  technique,  they  obtained  amplitude  variations  as  high  as  a 
factor  of  2.5.  Thus,  this  type  of  model  shows  some  promise  of  explaining  the  amplitude 
anomalies  from  NTS  explosions.  However,  it  would  be  preferable  to  invert  the 
amplitudes  themselves  for  a  candidate  structure,  rather  than  using  ad  hoc  models. 

Our  procedure  for  imaging  the  velocity  structure  uses  the  thin  lens  projection  to 
construct  a  model  of  a  warped  velocity  discontinuity.  Travel  time  delays  are  computed 
from  the  equation  of  wavefront  curvature  as  in  Haddon  and  Husebye  (1978),  and  are 
then  mapped  into  topography  of  a  velocity  discontinuity.  This  structure  is  then  tested 
by  computing  synthetic  seismograms  using  the  Kirchhoff-Helmholtz  integral  method. 

OBSERVED  TRAVEL  TIME  AND  AMPLITUDE  VARIATIONS 

Over  2000  P  wave  ab  (first-peak-to-first-trough)  amplitudes  and  travel  times  were 
measured  from  waveforms  for  57  explosions  at  the  Pahute  Mesa  and  Yucca  Flat  test  sites 
in  Nevada  (Fig.  1).  The  receivers  were  short-period  WWSSN  and  Canadian  Seismic 
Network  seismographs  at  distances  of  25*-95*.  A  complete  statistical  analysis  of  this 
data  set  is  presented  in  Lynnes  and  Lay  (1987a).  The  amplitude  data  set  actually 
consists  of  Amb  values  (i.e.,  deviations  from  the  explosion  mean).  The  station  average 
values  for  each  test  site  (pMSj  and  y^Sj)  are  presented  on  an  equal-area  diagram  in  Fig. 

2.  The  two  test  sites  have  similar  patterns  for  both  measurements,  suggesting  that  both 
the  amplitude  and  travel  time  patterns  have  a  component  related  to  deep  path  and/or 
near-receiver  effects.  However,  the  deviations  from  the  mean  for  Pahute  Mesa  are 
systematically  and  coherently  larger  in  each  case.  These  intersite  differences  indicate  a 
near-source  component.  In  an  attempt  to  minimize  the  deep  path  and  near-receiver 
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effects,  the  mean  of  the  average  station  anomalies  ([pMSi+YFSi]/2*NTSSi)  for  the  two  test 
sites  was  removed  from  the  individual  travel  time  and  amplitude  measurements.  This 
should  provide  a  fairly  conservative  estimate  of  the  near-source  component  of  amplitude 
and  travel  time  variation,  although  there  is  some  potential  for  effects  unique  to  one 

subsite  to  be  folded  into  the  data  for  the  other.  An  independent  test  of  this  assumption 

\ 

can  be  provided  in  the  case  of  the  amplitude  data  by  comparing  the  root-mean-square 
of  the  first  five  seconds  (RMS0*8),  which  represents  the  direct  arrival,  to  the  root-mean- 
square  of  the  next  ten  seconds  (RMS8*18),  which  represents  the  early  coda.  The  early 
coda  should  average  out  much  of  the  near-source  focussing-defocussing  effects,  thus 
providing  a  reference  for  the  direct  arrival  amplitudes.  The  station  mean  values  of 
log(RMS°*5/RMS8*18)  are  plotted  in  Fig.  3,  along  with  the  individual  station  mean 
magnitude  anomalies  with  the  regional  means  removed  (PMSi-NTgSl  and  ypS.-j^gS;).  The 
patterns  for  Pahute  Mesa  are  strikingly  similar  in  range  and  sense.  The  resemblance  of 
the  Yucca  Flat  patterns  is  somewhat  less,  probably  due  to  the  complex  very-near-source 
structure  (Lynnes  and  Lay,  1987b). 

Another  important  correlation  is  that  between  travel  times  and  amplitudes.  Fast 
travel  times  (negative  residuals)  generally  correspond  to  low  amplitudes,  e.  g.,  in  the 
north-northeast  of  the  Pahute  Mesa  pattern.  This  type  of  correlation  is  generally  an 
indication  of  focussing-defocussing  effects  (e.  g„  Chang  and  von  Seggern,  1980;  Butler, 
1983).  Thermally  activated  variations  in  attenuation,  on  the  other  hand,  would  predict  a 
correlation  of  slow  times  with  low  amplitudes  (cf.,  Butler,  1984). 

BACK  PROJECTIONS  OF  AMPLITUDES  AND  TRAVEL  TIMES 

The  procedure  for  imaging  thin  lens  structures  that  might  produce  the  observed 


amplitude  and  travel  time  patterns  is  a  simple  back- projection  of  the  teleseismic 
anomalies.  The  rays  are  traced  through  a  homogeneous  structure,  and  the  appropriate 
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anomaly  is  projected  to  the  point  below  the  source  area  where  the  ray  reaches  a  given 
depth  (Fig.  4a),  in  this  case  100  km.  The  structure  used  is  the  reference  model  for  NTS 
used  by  Taylor  (1983)  for  three-dimensional  travel  time  inversion,  on  top  of  the  T7 
model  of  Burdick  and  Helmberger  (1978).  The  projection  surface  is  then  gridded,  and 
the  values  within  each  grid  square  are  averaged  to  obtain  a  bin-averaged  projection  (Fig. 
4b).  The  bin-averaged  projections  for  depths  of  35,  100  and  160  km  are  shown  in  Fig. 

5  for  both  the  amplitude  and  travel  time  data  sets.  For  the  shallow  depth,  the  images 
formed  by  the  travel  time  and  amplitude  data  are  less  coherent  than  the  other  depths, 
with  mingled  negative  and  positive  values.  This  is  due  to  the  strong  similarity  of 
patterns  for  the  explosions  in  a  given  test  site:  a  structure  as  shallow  as  35  km  does  not 
affect  all  the  explosions  similarly,  leading  to  destructive  interference  in  the  back- 
projections.  At  100  km  depth,  the  images  become  more  coherent,  as  indicated  by  the 
larger  groupings  of  negative  values,  particularly  in  the  NNE.  Both  the  amplitudes  and 
travel  times  form  a  northeast-trending  low  under  the  test  site.  The  160  km  projection  is 
fairly  similar  to  the  100  km  projection 

A  similar  procedure  can  be  applied  to  the  deviations  from  the  subsite  station  mean. 
In  this  case,  both  the  event  mean  and  individual  station  mean  for  a  given  test  site  are 
subtracted  from  each  measurement.  At  large  depths,  the  anomalies  for  a  given  station 
will  fall  in  a  single  bin  and  average  out  to  zero.  Consequently,  this  component  of  the 
data  set  is  most  sensitive  to  very  shallow  structure. 

The  back-projections  of  the  individual  deviations  from  the  station  mean  are  shown 
for  Pahute  Mesa  in  Fig.  6.  The  shallowest  depth  has  a  relatively  incoherent  pattern, 
indicating  that  no  simple  structure  at  this  depth  can  explain  much  of  the  magnitude 
deviations.  However,  the  25  km  deep  projection  has  some  very  coherent  images, 
particularly  the  large  area  of  low  amplitudes  under  the  north-central  part  of  the  mesa. 
The  50  km  deep  projection  is  less  coherent  and  has  more  bins  with  zero  values,  as 
expected.  The  corresponding  projections  for  the  Yucca  Flat  test  site  do  not  in  general 
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form  coherent  images.  This  result  is  in  accord  with  the  hypothesis  of  Lynnes  and  Lay 
(1987b)  that  the  direct  arrivals  from  Yucca  Flat  are  heavily  influenced  by  very  shallow 
structure. 

CALCULATION  OF  TIME  DELAYS  FROM  AMPLITUDE  VARIATIONS 

The  criterion  we  use  for  choosing  a  projection  depth  is  the  correspondence  of  the 
travel  time  projection  with  the  travel  times  predicted  by  the  amplitude  projection, 
calculated  using  the  wavefront  curvature  equation  of  Haddon  and  Husebye  (1978).  The 
advance  (r(x,y))  of  a  plane  wavefront  at  an  infinitely  thin  lens  is  given  as 

-  hIr„r„-(rw)1l,  (1) 

where  A/Ao  -  relative  amplitude, 
h  -  depth  of  the  lens, 

and  subscripts  indicate  partial  derivatives  in  cartesian  coordinates. 

This  equation  is  essentially  a  non-linear  version  of  Poisson’s  equation  and  can  be  solved 
for  r(x,y)  given  sufficiently  smooth  boundary  conditions  and  amplitude  variations  by 
calculating  the  non-linear  terms  on  the  right  side  of  the  equation  from  the  previous 
iteration. 

Haddon  and  Husebye  (1978)  used  observed  travel  time  residuals  as  boundary 
conditions  in  their  study  on  the  NORSAR  array.  However,  for  an  irregular  geometry 
such  as  the  NTS  explosion  array,  this  choice  of  boundary  conditions  risks  introducing 
biases  based  on  a  small  subset  of  the  travel  times.  Furthermore,  this  subset  is  the  most 
poorly  determined  since  its  position  on  the  fringe  of  the  projection  leads  to  a  relatively 
low  number  of  rays  per  bin  and  precludes  a  significant  number  of  crossing  rays. 
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Consequently,  we  use  zero  boundary  conditions,  which  are  implemented  by  expanding 
both  the  travel  time  and  amplitude  grids  by  a  few  grid  points  and  smoothly  tapering 
each  to  zero  at  the  grid  edges. 

The  solution  to  equation  (1)  is  obtained  through  a  numerical  solution  to  Poisson’s 
equation  (Sweet,  1974;  Swartztrauber  and  Sweet,  1975)  on  a  rectangular  grid.  In  order  to 
satisfy  this  constraint,  data-poor  rows  and  columns  were  trimmed  from  the  edges  of  the 
projections  and  the  grid  was  expanded  2-4  points  to  implement  the  zero  boundary 
conditions  (Fig.  7a).  Values  were  then  interpolated  at  the  grid  points  corresponding  to 
bins  with  no  data  (Fig.  7b). 

A  spatial  filter  was  applied  to  both  the  amplitudes  and  travel  times.  The  grids  were 
transformed  using  a  2-dimensional  FFT,  and  the  higher  wavenumber  components  were 
removed  (Fig.  7c).  The  travel  time  surface  was  heavily  filtered  for  two  reasons.  Firstly, 
the  short- wavelength  components  are  less  reliable  than  the  long-wavelength  components 
due  to  station  effects  and  measurement  uncertainty,  and  are  less  stable  with  respect  to 
the  bin-averaging  procedure.  Secondly,  the  double  integration  of  the  amplitudes  in 
solving  equation  (1)  produces  perforce  a  very  smooth  travel  time  surface,  meaning  that 
the  amplitudes  cannot  resolve  short-wavelength  features.  Consequently,  the  amplitudes 
were  only  slightly  filtered  to  remove  short-wavelength  noise  before  solving  the  equation. 
The  spatial  filter  allows  a  comparison  between  features  of  similar  scale  length  between 
the  travel  time  and  amplitude-derived  time  delay  surfaces.  The  comparison  criterion  is 
a  cross-correlation  coefficient  of  the  filtered  travel  time  model  with  the  amplitude- 
derived  time  delay  model  at  the  grid  points  whose  corresponding  bins  contained  data. 

Models  derived  from  amplitude  variations 

The  results  for  the  thin  lens  models  which  had  the  best  match  of  amplitude- 
predicted  time  delays  and  observed  time  delays  are  shown  in  Fig.  8a.  The  depth  for 
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these  models  is  160  km,  and  a  correlation  coefficient  of  0.48  was  obtained.  The  models 
appear  most  similar  in  the  area  of  fast  travel  times  trending  northeast.  This  is  also  the 
area  of  the  projection  surface  with  the  most  data,  and  consequently  the  most  robust 
features.  If  the  observed  and  predicted  travel  time  surfaces  are  filtered  more  in  order  to 
compare  the  long-wavelength  features,  the  correlation  coefficient  increases  to  0.79  (Fig. 
8b).  The  amplitude  of  the  variation  in  the  two  models  is  also  comparable.  The  high 
correlation  in  the  long- wavelength  components  of  predicted  and  observed  travel  times  is 
a  strong  indication  that  most  of  the  robust  amplitude  variation  arises  from  elastic 
focussing-defocussing.  The  correlation  coefficient  peaks  at  160  km,  but  decreases  slowly 
with  increasing  depth.  This  is  because  the  back-projections  change  slowly  as  the  source 
array  appears  more  like  a  point  source.  However,  at  depths  below  160  km,  the 
amplitudes  predict  a  much  larger  range  in  travel  time  variation  than  indicated  by  the 
observed  time  delay  surface.  Therefore,  we  choose  160  km  as  the  optimal  projection 
depth. 

The  same  procedure  can  be  applied  to  the  results  of  the  deviations  from  the  subsite 
station  means.  For  Pahute  Mesa,  the  best  correlation  of  predicted  to  observed  travel 
times  is  at  a  depth  of  25  km,  which  also  had  the  most  coherent  images  for  the  amplitude 
scatter  (Fig.  6).  These  correlate  at  a  level  of  0.43,  while  the  very-long-wavelength 
component  correlates  at  a  0.71  level  (Fig.  9).  However,  the  predicted  range  in  travel 
times  is  only  0.05  s  for  the  very-long-wavelength  component,  whereas  the  observed 
range  is  0.20  s.  This  probably  indicates  that  a  single  thin  lens  cannot  predict  both  the 
amplitudes  and  travel  times  in  this  case,  but  rather  a  continuous  distribution  of  velocity 
heterogeneity  is  needed  to  accumulate  the  appropriate  travel  time  anomalies. 

Nevertheless,  it  is  instructive  to  see  how  far  a  discontinuity  derived  from  the  lens  can  go 
in  producing  the  observed  trends  in  amplitudes  across  the  mesa.  The  same  procedure 
was  attempted  for  Yucca  Flat,  but  the  predicted  travel  time  variations  were  vanishingly 
small  due  to  the  failure  of  the  back-projections  to  form  coherent  images. 
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The  three-dimensionally  varying  velocity  model  is  derived  by  simply  mapping  the 
time  delay  surface  into  vertical  deflections  of  a  velocity  discontinuity,  using  the  equation 

z(x,y)  -  At  /  (1/^  -  l/al),  (2) 

where  and  ctj  are  the  P  wave  velocities  above  and  below  the  velocity  discontinuity, 
respectively.  The  velocities  used  for  the  160  km  deep  case  are  7.70  km/s  and  8.25  km/s, 
chosen  to  correspond  to  the  range  of  velocities  in  the  low  velocity  zone  of  the  T7  model. 
We  model  the  shallow  (25  km)  case  as  a  Conrad- type  discontinuity  with  velocities  of  6.4 
km/s  and  7  km/s.  The  resultant  surfaces  are  shown  in  Fig.  20. 

Both  surfaces  are  characterized  by  broad  central  highs,  engendered  by  the  low 
amplitudes  in  the  center  of  the  projections,  with  deep  flanking  lows.  The  160  km  deep 
surface  has  a  range  of  about  30  km,  while  the  25  km  deep  surface  has  a  range  of  8  km. 
The  discontinuities  hypothesized  for  the  purposes  of  the  model  construction  are  the  base 
of  the  low-velocity  zone  and  the  base  of  the  upper  crust.  As  it  is  quite  plausible  that 
such  boundaries  would  be  rather  irregular,  the  topography  derived  herein  seems 
physically  reasonable. 

CALCULATION  OF  KIRCHHOFF- HELMHOLTZ  SYNTHETICS 

Haddon  and  Buchen  (1981)  and  Scott  and  Helmberger  (1985)  have  presented 
developments  of  the  three-dimensional  Kirchhoff-Helmholtz  integral  formulation  for  the 
transmitted  case.  The  equation  for  a  transmitted  potential  across  a  velocity  discontinuity 
E  is: 

<P2  •  l/(4ff)/aE  T  (’(t-T/aj+rJaJ  (cos©2  +  cos01)/(rora2)  dS, 
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where  T  is  the  acoustic  plane  wave  transmission  coefficient, 
r  is  the  source-surface  element  distance, 

O 

r  is  the  surface  element- receiver  distance, 

is  the  angle  between  the  surface  and  the  refracted  wavefront, 

0j  i3  the  angle  between  the  surface-receiver  vector  and  the  normal  to  the  surface, 
and  dS  is  an  areal  element  (Fig.  11). 

(This  is  modified  from  equation  13  in  Scott  and  Helmberger,  in  which  the  second  cosine 
term  should  be  divided  by  aJt  rather  than  a1.)  The  transmission  coefficient  is  calculated 
for  the  local  angle  of  incidence  on  the  surface.  If  this  angle  exceeds  the  critical  angle, 
the  transmission  coefficient  is  set  equal  to  zero.  A  synthetic  ramp  response  is  calculated 
by  simply  summing  step  functions  with  the  appropriate  multipliers  and  time  lags  for 
each  surface  element.  In  order  to  obtain  a  teleseisraic  displacement  seismogram,  the 
resulting  ramp  response  is  numerically  differentiated  (Fig.  12),  and  convolved  with  the 
second  derivative  of  an  explosion  source  potential,  a  Futterman  operator  with  t*  of  0.75, 
and  a  short-period  WWSSN  instrument  response.  The  source  used  for  these  seismograms 
is  a  modified  Haskell  source  (Helmberger  and  Hadley,  1981)  with  b-1  and  K=8.  A  pP 
phase  with  a  time  lag  of  0.9  s  and  reflection  coefficient  of  -1  is  also  added,  but  does  not 
affect  the  ab  part  of  the  waveform. 

The  seismograms  are  calculated  for  two  welded  half-spaces,  using  acoustic 
transmission  coefficients.  The  receivers  are  located  20,000  km  below  the  surface, 
corresponding  to  the  geometric  spreading  at  about  60°  in  the  real  earth,  and  at  the 
appropriate  x-y  coordinates  so  that  the  intersection  of  the  geometric  ray  with  the  flat 
surface  coincides  with  the  appropriate  site  in  the  back-projections.  Synthetic 
seismograms  are  also  computed  for  the  flat-layer  case.  The  ab  amplitudes  are  calculated 
for  the  model  synthetics  and  divided  by  the  ab  amplitudes  for  the  flat-layer  synthetics. 

In  all,  the  computation  of  the  synthetics  is  fairly  easy  and  quick,  facilitating  the  testing 
of  a  large  number  of  different  models  for  many  source-station  pairs. 
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Synthetic  tests  of  models 

For  the  160  km  deep  discontinuity,  synthetic  seismograms  were  calculated  for  56 
teleseismic  stations  from  sources  in  the  center  of  the  Pahute  Mesa  and  Yucca  Flat  sites. 
The  differences  in  the  Pahute  Mesa  and  Yucca  Flat  log  amplitudes  are  shown  for  each 
station  in  Fig.  13,  along  with  the  differences  in  the  observed  Pahute  Mesa  and  Yucca 
Flat  station  mean  magnitudes  (p^-y^S;).  The  patterns  are  strikingly  similar,  differing 
only  in  strength.  Even  fairly  subtle  details  in  the  observed  pattern  are  matched  in  the 
synthetics,  such  as  the  trend  in  the  closer  stations  from  relatively  low  magnitudes  to  the 
north  for  Pahute  Mesa  relative  to  Yucca  Flat,  to  high  magnitudes  in  the  NNE,  low 
magnitudes  in  the  NE,  and  high  magnitudes  in  the  ENE  and  east.  Thus,  the  imaging  of 
three-dimensionally  varying  structure  from  the  amplitude  variations  appears  to  be 
successful  in  predicting  much  of  the  amplitude  behavior.  Of  course,  it  cannot  be 
expected  to  match  the  full  range  of  variation  because  the  back-projection  utilizes 
destructive  interference  to  mask  out  anomalies  that  are  not  identified  with  that  particular 
projection  depth.  Thus,  the  technique  is  intrinsically  conservative. 

The  patterns  of  travel  time  residuals  also  are  strikingly  similar  for  both  the  synthetic 
and  observed  seismograms  (Fig.  14).  Again,  the  range  of  variation  is  not  matched  by  the 
model,  due  to  the  conservative  nature  of  the  procedure.  Nevertheless,  it  is  encouraging 
that  the  model  derived  from  amplitude  variations  qualitatively  predicts  the  observed 
travel  time  variations  as  well.  This  agreement  also  indicates  that  the  choice  of  zero 
boundary  conditions  in  the  wavefront  curvature  equation  does  not  introduce  any 
significant  artifacts. 

The  Kirchhoff-Helmholtz  synthetics  for  the  shallow  model  are  somewhat  more 
difficult  to  construct.  Because  of  the  shallow  depth  and  short  scale  lengths  of  the 
model,  the  local  angle  of  incidence  exceeds  the  critical  angle  for  part  of  the  surface  (Fig. 
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15a)  for  many  of  the  source-station  geometries.  In  fact,  this  effect  is  a  primary 
contributor  to  the  variations  in  waveform  and  amplitude  for  the  different  structures 
tested  by  Scott  and  Helmberger  (1985).  For  the  model  tested  here,  the  masking  effect  is 
particularly  prevalent  on  the  steep  transition  from  the  high  to  the  low  on  the  surface. 
Since  such  surface  elements  are  assigned  a  transmission  coefficient  of  zero  in  the 
computation  of  the  seismograms,  this  results  in  a  somewhat  artificial  truncation  effect  in 
the  very  early  part  of  the  waveform  (Fig.  15b).  These  truncation  effects  can  produce 
large  variations  in  the  ab  amplitudes.  In  order  to  assess  the  general  magnitude  of  these 
effects,  the  synthetics  were  recomputed  by  assigning  supercritical  elements  the  amplitude 
factor  of  adjacent  elements.  This  has  the  effect  of  generally  smoothing  out  the  effect  of 
the  truncation  (Fig.  15c).  Stations  for  which  the  ab  amplitudes  changed  significantly, 
which  were  primarily  the  stations  to  the  west  such  as  MAT,  were  excluded  from  the 
analysis.  However,  the  actual  ab  amplitudes  used  herein  are  those  calculated  by  masking 
out  critical  and  supercritical  elements. 

The  variation  of  log  amplitude  with  position  in  the  test  site  for  three  stations  are 
shown  in  Fig.  16.  Since  the  corresponding  variations  in  observed  mb  values  are  highly 
contaminated  by  incoherent  measurement  errors,  the  deviations  from  the  station  means 
for  nearby  stations  were  averaged  in  for  each  explosion  for  the  comparison  (Fig.  16).  As 
with  the  deep  model,  the  synthetic  log  amplitudes  have  the  correct  sense  of  variation, 
but  underpredict  the  amplitudes.  However,  this  qualitative  correlation  is  quite 
important,  because  it  implies  that  the  "scatter"  about  the  statior  means  can  in  fact  be 
partially  interpreted  in  a  deterministic  framework.  This  is  supported  by  the  coherent 
variations  of  observed  mb  deviations  within  the  mesa  for  each  group  of  stations. 


DISCUSSION 
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It  is  interesting  to  compare  the  models  derived  from  amplitudes  in  this  study  with 
those  derived  from  travel  times.  A  vertical  cross-section  through  the  Minster  et  al. 
(1981)  and  Taylor  (1983)  structures  of  NTS  is  shown  in  Fig.  17,  with  cross-sections  of 
the  two  velocity  discontinuity  models  derived  from  amplitudes  superimposed.  The  two 
travel  time  structures  are  reasonably  similar,  although  the  Taylor  (1983)  structure 
terminates  more  shallowly  and  has  much  larger  velocity  variations.  The  amplitude- 
derived  structure  at  25  km  is  qualitatively  similar  to  both,  with  an  upwarp  in  the  high- 
velocity  material  closely  corresponding  to  fast  velocity  regions  in  both  travel  time 
models.  The  deep  amplitude  model  is  more  difficult  to  compare  since  it  is  located 
deeper  than  either  of  the  travel  time  model  parameterizations.  However,  if  the  high- 
velocity  regions  in  the  travel  time  models  are  extended  down  along  the  raypaths  to  the 
stations  to  the  NNE,  they  would  intersect  an  upwarp  of  the  high-velocity  material. 

The  comparison  is  complicated  by  the  dichotomy  of  velocity  boundary  undulation 
models,  which  concentrate  velocity  heterogeneity  over  a  small  depth  range,  and 
continuous  block  models,  which  have  a  tendency  to  smear  the  heterogeneity  along  the 
dominant  raypaths.  The  relative  ability  of  each  class  of  model  to  produce  amplitude 
variations  for  a  given  set  of  travel  time  perturbations  is  contingent  on  whether  log 
amplitudes  accumulate  linearly  along  the  raypath.  Haddon  and  Husebye  (1978) 
computed  amplitude  anomalies  for  a  150  km  deep  thin  lens,  and  a  100  km  thick  lens 
models  under  NORSAR  that  produces  the  same  propagation  times  for  vertical  incidence. 
The  patterns  were  nearly  identical,  suggesting  log-linear  accumulation  of  amplitude 
anomalies.  However,  as  noted  above,  the  model  of  Taylor  (1983)  for  NTS  produces  much 
larger  variations  that  the  152  km  thick  model  of  Minster  et  al.  (1981),  despite  the  fact 
that  they  satisfy  similar  travel  time  residual  data  sets.  Thus,  it  seems  that  a  log-linear 
accumulation  model  for  amplitudes  may  be  valid  at  depth,  but  may  break  down  near  the 
source.  Indeed,  this  may  explain  why  the  thin  lens  amplitude  model  at  25  km  depth 
grossly  underpredicted  the  range  of  travel  time  variation.  Specifically,  if  the 
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heterogeneity  is  concentrated  over  a  small  depth  range  near  the  source,  a  much  smaller 
heterogeneity  is  necessary  to  produce  a  given  amplitude  range. 

One  shortcoming  of  the  procedure  used  in  this  study  for  deriving  three- 
dimensionally  varying  models  from  amplitudes  appears  be  a  systematic  underprediction 
of  the  range  in  amplitude  variations.  This  is  a  common  problem  in  attempts  to  predict 
amplitudes  in  an  elastic  framework.  As  noted  earlier,  the  structures  derived  from  travel 
times  by  Thomson  and  Gubbins  (1982)  for  NORSAR  and  by  Minster  et  a!.  (1981)  for 
NTS  both  underpredict  amplitude  variations  to  a  varying  degree.  The  structure  derived 
by  Taylor  (1983)  for  NTS  does  predict  amplitude  variations  of  the  same  order  as  the 
data  (Cormier,  1987),  but  this  is  the  result  of  very  large  velocity  differences  over  short 
distances,  ranging  from  16%  over  14  km  in  the  shallowest  layer  to  5%  over  20  km  in  the 
deepest  layer. 

In  theory,  a  structure  derived  from  the  amplitudes  should  be  better  able  to  predict 
the  range  of  amplitude  variation.  Thomson’s  (1983)  models  for  NORSAR  derived  using 
a  ray- theoretical  amplitude  inversion  predict  the  observed  amplitude  variation  better 
than  the  travel  time  model  A2  (Thomson  and  Gubbins,  1982),  even  though  the  latter 
model  has  stronger  velocity  variations.  For  the  case  of  the  procedure  presented  in  this 
paper,  the  main  factor  contributing  to  the  underprediction  of  the  amplitudes  seems  to  be 
the  parameterization  of  only  a  single  surface  for  a  given  model.  The  destructive 
interference  in  the  back-projections  ensures  that  only  part  of  the  available  information 
in  the  amplitudes  is  being  utilized.  We  have  attempted  to  circumvent  this  by  examining 
both  the  regionally  demeaned  anomalies  (for  the  deep  models)  and  the  deviations  from 
the  individual  station  means  (for  the  shallow  models).  However,  it  may  be  possible  to 
simultaneously  back-project  anomalies  to  two  or  more  surfaces,  i.e.,  tomographically 
image  the  amplitude  data.  This  should  retain  more  of  the  amplitude  information, 
resulting  in  more  accurate  models. 
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It  should  also  be  noted  that  we  have  attempted  to  explain  only  the  intrasite  and 
intersite  differences.  However,  it  is  apparent  from  the  similarity  of  the  Yucca  Flat  and 
Pahute  Mesa  azimuthal  patterns  that  the  two  subsites  share  a  large  common  component 
in  both  amplitude  and  travel  time  variations  (Fig.  2),  with  low  magnitudes  and  fast 
travel  times  to  the  northeast.  This  component  is  as  much  as  0.2  units  for  the  magnitudes 
and  1.0  s  for  the  travel  times.  The  relatively  small  ratio  of  amplitude  to  travel  time 
anomaly  indicates  very  broad  scale  lengths,  and  this,  coupled  with  the  fact  that  both 
subsites  must  be  affected  similarly,  indicates  a  deep  (>400  km)  arce  for  the  common 
patterns. 

CONCLUSIONS 

We  have  used  a  combination  of  thin  lens  projections,  calculation  of  time  delays  from 
amplitudes  using  the  wavefront  curvature  equation,  and  mapping  time  delays  into  a 
warped  velocity  discontinuity  to  image  three-dimensionally  varying  structure  from 
amplitude  data.  Synthetic  seismograms  can  be  easily  computed  for  these  discontinuities 
using  the  Kirchhoff-Helmholtz  technique.  A  comparison  of  observed  and  synthetic 
magnitude  anomalies  indicates  that  the  coherent  intersite  and  intrasite  differences  in 
magnitudes  from  NTS  explosions  can  be  partially  explained  by  warped  velocity 
discontinuities  25  km  and  160  km  beneath  the  source  array. 
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Fig.  1.  Map  of  the  Nevada 
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Fig.  3.  Equal-area  plots  of  station  mean  magnitude  anomalies,  with  NTS  station  means 
removed,  are  plotted  above  station  mean  values  of  Iog(RMS°'s/RMS5*15). 
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Fig.  4.  Example  of  back- projection  procedure.  Magnitude  anomalies,  with  NTS  station 
means  removed,  are  projected  through  laterally  homogeneous  structure  to  a  depth  of  100 
km  (a).  The  surface  is  gridded  and  the  values  in  each  grid  square  are  averaged  to 
provide  a  bin-averaged  projection  (b).  The  triangles  indicate  the  positions  of  the 
sources. 
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Fig.  6.  Bin-averaged  back-projections  of  deviations  from  individual  station  means  for 
the  Pahute  Mesa  test  site  at  depths  of  10  km,  25  km  and  50  km.  A  coherent  image  of 
low  values  forms  in  both  the  magnitude  and  travel  time  projections  at  25  km  depth. 
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Fig.  7.  Example  of  processing  of  bin-averaged  surface  to  obtain  a  continuous,  filtered, 
rectangular  grid,  for  a  depth  of  160  km.  The  original  grids  are  shown  in  Fig.  5 
(bottom).  The  grid  is  cropped  to  remove  data- poor  columns  and  rows  and  expanded  to 
allow  the  implementation  of  zero  boundary  conditions  (a),  interpolated  to  obtain  a 
continuous  rectangular  grid  (b),  and  filtered  in  doubly -transformed  domain  (c).  The 
travel  times  are  more  heavily  filtered  than  the  amplitudes  (see  text).  Only  points  which 
have  data  are  shown  in  the  last  plot 


Fig.  8.  Thin  lens  results  for  160  km  depth,  filtered  to  retain  long- wavelength  (a)  and 
very-long-wavelength  (b)  components.  The  top  plot  plots  are  the  filtered,  observed 
travel  times,  and  the  bottom  are  the  filtered  travel  times  predicted  from  the  magnitude 
variations. 
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Fig.  10.  Contour  plou  (in  km)  of  the  velocity  discontinuity  surfaces  obtained  from 
inversion  of  amplitudes.  The  deep  (160  km)  model  obtained  from  intersite  differences 
in  magnitude  patterns  is  shown  on  the  left.  The  shallow  (25  km)  surface  obtained  from 
intrasite  differences  is  shown  on  the  right. 
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Fig.  12.  Example  of  synthetic  seismogram  calculation.  On  the  left  is  a  step 
(differentiated  ramp)  response  for  Pahute  Mesa  to  UME  for  the  deep  model.  This  is 
convolved  with  a  function  containing  the  source,  instrument,  free-surface  and 
attenuation  (t*— 0.75)  responses  (center).  On  the  right  is  the  final  seismogram  with  the 
ab  amplitude  relative  to  the  flat  layer  case  indicated  in  the  upper  right. 
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Fig.  13.  Comparison  of  synthetic  log  amplitude  variations  with  observed  variations  in 
mb.  The  equal-area  plot  on  the  left  shows  the  observed  station  average  mb  values  for 
Yucca  Flat  differenced  from  the  corresponding  station  means  for  Pahute  Mesa.  The 
equal-area  plot  on  the  right  shows  synthetic  log  amplitude  variations  of  Yucca  Rat 
explosions  differenced  from  those  of  Pahute  Mesa. 


Fig.  15.  The  25  km  deep  model  is  shown  in  (a),  the  stippled  area  indicating  where  the 
local  angle  of  incidence  exceeds  the  critical  angle  for  the  event  CAMEMBERT  (filled 
triangle)  to  station  MAT.  The  geometric  arrival  intersects  the  surface  at  the  star.  The 
'masked'  portion  of  the  surface  leads  to  a  truncation  in  the  first  second  of  the  step 
response  (b).  This  is  smoothed  out  by  assigning  the  masked  elements  the  amplitude 
factors  of  adjacent  elements  (c)  in  order  to  evaluate  the  magnitude  of  this  effect.  The 
numbers  in  (b)  and  (c)  indicate  ab  amplitudes  relative  to  the  flat-layer  case. 
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Fig.  16.  (a)  Variations  in  synthetic  log  amplitudes  with  position  in  Pahute  Mesa, 
calculated  for  Resolute  Bay,  Canada  (RES),  Umea,  Sweden  (UME),  and  Porto,  Portugal 
(PTO).  (b)  Observed  variations  in  mb,  obtained  by  averaging  nearby  stations  with  those 
for  which  the  synthetics  were  calculated. 
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Fig.  17.  Comparison  of  discontinuity  models  with  block  inversion  models  derived  from 
travel  times  by  other  authors.  The  figure  on  the  left  is  a  modification  of  Fig.  4.13b  of 
Minster  et  al.  (1981),  which  is  a  vertical  cross-section  through  their  model,  represented 
by  the  velocity  contours.  Symbols  indicate  the  velocity  perturbations  in  the  blocks  of 
Taylor’s  (1983)  model  that  the  cross-section  intersects.  The  thick  black  lines  are  the 
equivalent  cross-section  through  the  shallow  (25  km)  velocity  discontinuity  model  and  an 
extended  cross-section  through  the  deep  (160  km)  model.  The  locations  of  the  cross- 
sections  are  nlotted  on  contour  maps  of  the  surfaces  on  the  right.  High  velocity  regions 
or  extensions  thereof  in  the  block  travel  time  models  correspond  to  upwarps  of  high 
velocity  material  in  the  amplitude-derived  discontinuity  models. 
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